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1. INTRODUCTION 


Although from a historical point of view, composite materials have found practical use 
for centuries, during the past two decades there has been a tremendous increase in the use 
of composite materials in engineering applications, particularly in aerospace engineering. 
The main attraction of these materials is the ability to design and manufacture such 
materials to sustain a specific type of loading in the most efficient manner. If properly 
produced, composites can often achieve a combination of properties that are far superior 
to the properties of the individual constituents acting independently. 

It is well known that gas turbine engine structures, particularly those components di- 
rectly in the hot gas flow path, are subjected to extremely severe thermal and mechanical 
loading that can often lead to creep enhanced distortion, cracking and low cycle fatigue. 
As the demand for more efficient propulsion systems rise so does the thermal and mechan- 
ical loading. It is unlikely that the current generations of metal alloys would be suitable 
candidates for structural components in the future generations of efficient propulsion sys- 
tems. Ceramic components are often thought to be ideal as far as their thermal durability 
is concerned. Unfortunately, ceramics do not have adequate tensile strength to sustain a 
high level of mechanical loading. In recent years there has been significant effort in the 
attempt to incorporate fibrous inclusions within a ceramic matrix to develop a class of new 
materials, known as ceramic composites, for advanced engineering application. 

The mechanical behavior of ceramic composites under nonlinear, thermal, and dynamic 
loading is extremely complex and can only be understood if the observed behavior is 
interpreted in terms of micromechanical analyses. Such analyses must take care of the 
complex interaction of the individual fibers or bundles of fibers embedded in the three- 
dimensional ceramic matrix and must allow for increasing levels of sophistication in terms 
of the idealization of the fibers as well as the ceramic matrix. In addition complex interface 
behavior and controlled failure of the fiber must be considered. 
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It is evident that for proper micromechanical analysis of ceramic composites one needs 
to use a numerical method that is capable of idealizing the individual fibers or individual 
bundles of fibers embedded within a three-dimensional ceramic matrix. The analysis must 
be able to take account of high stress or temperature gradients from diffusion of stress or 
temperature from the fiber to the ceramic matrix and allow for the interaction between 
the fibers through the ceramic matrix. The analysis must be sophisticated enough to 
deal with failure of fibers described by a series of increasingly sophisticated constitutive 
models. Finally, the analysis must deal with micromechanical modeling of the composite 
under nonlinear thermal and dynamic loading. 

The boundary element method is uniquely suited for the task. BEM has proven its 
ability to accurately determine stress near a stress concentration. All functional quantities 
in a BEM system are on the boundary and interface surfaces, therefore, allowing nonlinear 
interaction on the interface between the matrix and the fiber to be readily described by fail- 
ure models. Furthermore, recent development has shown the generality and versatility of 
boundary element method in analyzing large two- and three-dimensional models subjected 
to static, dynamic, and thermal loads involving materials with nonlinear behavior. 

This report details progress made during the first three years of a five-year program 
commencing in March 1988, towards the development of a boundary element code designed 
for the micromechanical studies of advanced ceramic composite. Additional effort has been 
made in generalizing the implementation to allow the program to be applicable to real 
problems in the aerospace industry. 

The ceramic composite formulations developed for this work have been implemented in 
the three-dimensional boundary element computer code ‘BEST3D’ which was developed 
for NASA by Pratt and Whitney and the State University of New York at Buffalo under 
contract NAS3-23697. BEST3D was adopted as the base for the ceramic composite pro- 
gram, so that many of the enhanced features of this general purpose boundary element 
code could be utilized. Some of these facilities include sophisticated numerical integration, 
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the capability of local definition of boundary conditions, and the use of quadratic shape 
functions for modeling geometry and field variables on the boundary. The multi-region 
implementation permits a body to be modeled in substructural parts; thus dramatically 
reducing the cost of the analysis. Furthermore, it allows a body consisting of regions of 
different ceramic matrices and inserts to be studied. 

In the next chapter the elastostatic and steady-state heat conduction boundary ele- 
ment formulation for ceramic composites is developed. The method for the semi-analytic 
integration of the kernel functions about the fiber is presented and the numerical imple- 
mentation of the formulation is discussed. This is followed by the development of the 
uncoupled thermoelastic BEM formulation in Chapter 3. Up to this point, only fibers 
assumed to be perfectly bonded to the composite matrix were considered. In Chapter 4, 
the previous formulations are rederived in a form suitable for more general fiber to ma- 
trix interface connections. A spring and a nonlinear, coulomb friction constitutive models 
are developed for use as interface relations, and nonlinear solution algorithm is presented. 
Chapter 5 contains the transient heat conduction and transient uncoupled thermoelastic 
BEM formulations. Chapter 6 outlines the development of the general computer program 
implementation. In Chapter 7, a number of numerical examples Eire presented to demon- 
strate the power of the present implementations. This report is then concluded with a 
summary and plan for future developments. 
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2. Steady-State Heat Conduction and Elastostatic BEM Formulation 

2.1 Introduction 

The conventional boundary integral equations for elastostatic and steady-state heat 
conduction analyses are used in deriving a boundary element formulation for the analysis 
of ceramic composite structures. The boundary integral equation written for a point in the 
interior of the composite matrix is modified by adding to it the boundary integral equations 
of each insert written at the same point in the composite matrix. This eventually eliminates 
the displacement (or temperature) variables on the insert-matrix interface from the system, 
and therefore, reduces the total number of equations required for a solution of the system. 

The BEM formulation for the steady-state heat conduction analysis of ceramic com- 
posites is identical to the elastostatic formulation, and therefore, the two formulations will 
be derived as one. 


2.2 Boundary Integral Equation Formulations 

The direct boundary integral equation for the displacement (or temperature) at a point 
£ inside an elastic composite matrix is 


<?<*(0«<(0 = J [q§(*,0<?(*) - F£(*,0“?(*) dS{x) 

n=l Js * *- 


dS"(x) 


( 2 . 1 ) 


i,j = 1,2,3 for elastostatics 
i,j = 1 for heat conduction 


where 

Gij , Fij are the fundamental solutions of the governing differential equations of the ceramic 
matrix of infinite extent 

Cij axe constants determined by the geometry at £ 

U{,ti are displacements and tractions (or temperature and flux) 


4 


5, S" are the surfaces of the outer boundary of the matrix and the n th hole (left for fiber), 


respectively 

N is the number of individual insert fibers 


Superscripts O and H identify the quantities on the outer surface of the matrix and the 
quantities on the surface of the hole, respectively. 

The conventional boundary integral equation for displacement (or temperature) can 
also be written for each of the N insert fibers. For the displacement (or temperature) at a 
point £ inside the n th insert we can write 


0 = J sn 


dS n (x) 


( 2 . 2 ) 


i,j = 1,2,3 for elastostatics 
i, j = 1 for heat conduction 


G/j , Fjj are the fundamental solutions of the n th insert 

C(j are constants determined by the geometry at { in insert n 
u!,t! Eire displacement and tractions (or temperature and flux) associated with the n th 
insert 

5" the surface of the n th insert 


We next examine the interface conditions between the composite matrix and the insert. 
For a perfect bond the displacement (or temperature) of the matrix and the displacement 
(or temperature) of the inserts are equal and the tractions (or fluxes) along the interface 
are equal and opposite. 

uf (x) = u,-(x) (2.3a) 

if (x) = -tf(x) (2-3 b) 

In elastostatics, when the eleistic modulus of the insert is much greater than the modu- 
lus of the composite matrix, the Poisson ratio of the insert can be assumed equal to that of 
the matrix with little error (no approximation is required for the heat conduction analysis). 
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Therefore, upon consideration of the surface normals at the interface and examination of 
the Fij kernels, we can write the following relation for the n th insert 


Fij(x,Z) = -F%(x,t) 


(2.3c) 


Substitution of equations (2.3) into equation (2.2) yields the following modified boundary 
integral equation for insert n. 

C£(0m( 0 = J sn [- G lj( x ’Oii( x ) + (*.0«f(*)] ds "(*) ( 2 - 4 ) 

Finally adding the N insert equations (2.4) to equation (2.1) and cancelling terms, yields 
the modified boundary integral equation for the composite matrix 


0 = j s [og(*,0<f(x) - F^(*,o«f(*) 

+ E/ GZ(x,Ot?(x)dS"(x) 


dS(x) 


(2.5) 


where 

^ij( z >0 = G ij( x >0 ~ ( G ij( x >0) n 

Cij(f) axe constants dependent on the geometry for a point £ on the outer boundary and 
Cij({) = Sij for a point ( in the interior of the body. 


2.3 Analytic Integration Around an Insert 

The boundary element discretization of equations (2.4) and (2.5) in the conventional 
manner [Banerjee and Butterfield (1981)] requires a very fine discretization about the 
hole/insert. Alternatively, a new formulation is introduced in this report for the efficient 
modeling and analysis of holes/inserts using what the authors refer to as ‘Insert Elements’. 
The inserts are defined with Insert Elements by describing the centerline of the (curvilin- 
ear, tubular) insert with nodal points; defining the connectivity of the nodal points; and 
specifying the radius of the insert at each of these nodal points. Internally the program 
generates the surface of the insert and the hole in which the resulting displacements (or 
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temperatures) and tractions (or fluxes) are described using a trigonometric circular shape 
function in the circumferential direction and a curvilinear shape function of any order 
in the longitudinal direction (the present work employs both linear and quadratic shape 
functions for this purpose). A long hole (which is allowed to vary in diameter) can be 
described by a number of the insert elements connected end to end, and any insert element 
not connected to another is assumed, by the program, to be closed at the end by a circular 
disc. 

Using the concept of the insert element, the essential part of the formulation is the 
conversion of the two-dimensional surface integration of the insert (and of the hole) to a 
one-dimensional integration. By performing a semi-analytical integration on the surface 
of the hole (or insert) the numerical integration is significantly reduced. In equation (2.5) 
the integral under the summation is the integral associated with the hole which is to be 
modified. To facilitate an analytic integration in the circumferential direction, the three- 
dimensional kernel functions are first expressed in local coordinates with the center of the 
coordinate system coinciding with the center of the insert/hole and the z axis aligned with 
the centerline of the insert. The relative translation £,■ is added to the field coordinate & 
and the rotation is applied using the appropriate vector transformation. 

— a ij£j 

where aij are the direction cosines between the axis of the local and global coordinate 

systems and the bar indicates a local variable. 

The integration point Xi for a ring can now be expressed in cylindrical coordinates 

relative to the center of the hole/insert as 

x\ = RcosO 

X2 = RsinO 

X3 = 0 

where R represents the radius of the insert, i.e., R = {x\ + 
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The normal vectors are transformed by 


u\ = n r cos6 


n-2 = n r sinO 


Tir 3 — Tl z 

where n r and n z represents the normals of the side of the hole in local coordinates and are 
dependent on the change in the radius of the insert /hole. On the side of a straight hole 
n r - 1 and n z = 0, and on the flat surface closing the end of the hole/insert n r = 0 and n z — 1. 

Next a circular shape function is employed to approximate the variation in the trac- 
tion (or flux) about the circumference of the hole/insert. The circular shape function is 
multiplied and integrated with the three-dimensional kernel, allowing the nodal values 
of traction (or flux) to be brought outside the integral. The shape function is expressed as 

ti = M 7 *7 (summation over 7 is implied, 7 = 1, 2, 3 ) 

where 

M\&) = i + ^cos6 

M 2 (6) = ^ + —■ sinB — i cos$ 

1 v/S 1 

M 3 (6) = —sin8 — -cosd 

OO u 

and q is the nodal traction (or flux). 

A modified circular shape function is used in the integration over the end of the hole to 
insure continuity of traction at the center of the end surface. The modified shape function 
is expressed as: 

M 7 = aM 7 + fc/3 7 = 1.2,3 

a = r/R b = (R — r)/R 

where 

R is the radius of the hole at the end, 

r is the location of the integration (Gauss) point as it sweeps from 
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r = 0 to r = R, and 

M 7 is the circular shape function defined above. 

The traction must also be transformed between the local and the global systems by 

ti ~ dik^k 


or 

tj — d m jt m 

Since flux is a scalar quantity, no transformation is necessary in the heat conduction 
analysis. 

The last term in equation (2.5) can now be analytically integrated in the circumferential 
direction. For the m th hole the two integrals involved can be expressed as 

I GUx,Z)t?(x)ds m (x)= f a jk I*" G l £*\R,e,z,()M-iRdea it qdC m 
J s m Jc m Jo 

= f cr? j (R,z,oqdc m (z) 

where the indicated integration over C m is now a one- dimensional curvilinear integration 
along the hole and GR represent the analytically integrated hole/insert kernels. Note, since 
the transformation vector a ik is the independent of angle 0, it may be taken outside the dO 
integration. Similar analytic integration is also performed in equation (2.4). 

The final kernel functions of the hole/insert obtained from the analytical integration 
are very lengthy and therefore are not presented in this report. They contain functions of 
elliptical integrals which in general are expressed numerically by common series approxi- 
mations. For a range of input values (coordinate locations), several higher order elliptical 
integral functions were found to produce incorrect numerical results. To overcome this 
problem, several new series were derived using a best fit polynomial approximation (as a 
function of the modulus of elliptic integrals) using values of the integrals calculated by a 
very accurate numerical integration in the circumferential direction. 
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The derivation of the insert /hole kernels corresponding to the strain equation are 
accomplished from the displacement equation (2.5) by differentiation and application of 
the strain-displacement equations. The stress equation is then found using Hooke’s Law. 
Due to the complexity of the resulting insert/hole kernels, the authors perform the required 
differentiation before the analytic circumferential integration. 

Finally we note, an insert which has curvature along its length will differ in surface 
area about the circumference on the curved portion of the insert. This is neglected in the 
formulation since the analytical integration is performed on an axisymmetric ring in which 
the surface area is constant about the circumference. This error, however, is small and 
disappears completely on a straight tubular insert which is most commonly encountered. 
We should also mention that the inserts should not intersect the outer surface of the 
body or intersect other inserts. This minor restriction can be ignored if results at these 
locations are not of interest. As long as the inserts do not coincide with nodal points of 
other elements the errors will be localized and will not affect the overall boundary element 
solution. 


2.4 Numerical Implementation 

The integral representations of Section 2.2 are extremely accurate statements of the 
ceramic composite problem, however, approximations such as finite discretization and nu- 
merical integration are necessary in order to obtain a solution to non-trivial problems. 
The goal of the numerical implementation of the present formulation is to obtain the most 
accurate and efficient implementation possible. 

2.4.1 Discretization 

After the analytical integration in the circumferential direction is complete, an insert in 
a three-dimensional solid can be modeled as a two-dimensional curvilinear line element with 
a prescribed radius at each longitudinal node. In the present work, linear and quadratic 
shape functions are utilized in modeling the geometry and field variables along the insert 
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element as well as the boundary elements on the outer surface of the 3-D body. In the 
discretized form the displacement boundary integral equation for an elastic body containing 
inserts (equation (2.4)) can be expressed for a single insert as 


^(O^O — Ef/ GliMN^dCr 


i 7 


( 2 . 6 ) 


where 

P is the number of line elements, and 

represents a shape function over the curvilinear line element. Summation over y 


is implied. 

<7 and «7 are nodal values of traction (or flux) and displacement (or temperature) on 
the surface of the hole, respectively. 


In a similar manner, equation (2.5) can be discretized using one- and two-dimensional 
shape functions in the following manner. 

Q 


0 = £[ / Gf j (x,OL p (m,m)dsAt p i 

9=1 •"te' J 

~i\ f F t P(x,OL 0 (vum)dSO 

,= 1 L Js '' 

+ E[/ Gijix^WWdCr 
t^i 1 -Jc* 


( 2 . 7 ) 


where 


Q is the number of surface elements on the outer surface of the composite matrix in 
the region, and 


L 0 {i}\,r\ 2 ) represents a two-dimensional shape function. Summation over 7 and /? is 
implied. 


It is important to note that the displacement and traction (or temperature and flux) 
on an insert/hole varies in the longitudinal as well as the circumferential direction, i.e., for 
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displacement (or temperature). 

uj = M°N~ < u'p 

The circular shape function M a has been analytically integrated into the kernel func- 
tions of equations (2.6) and (2.7). The ends of the inserts are assumed to be a flat surface 
and a one- dimensional numerical integration is carried out in the radial direction. The 
coefficients obtained from the integration over the end are lumped with their respective 
coefficients from the integration of the side of the insert. 

Equations for stress, strain, or flux at points in the interior of a body can be discretized 
and integrated in a similar manner. 

2.4.2 Numerical Integration 

The complexity of the integral in the discretized equation necessitates the use of nu- 
merical integration for their evaluation. The steps in the integration process for a given 
element is outlined below: 

1. Using appropriate Jacobian transformations, a curvilinear line element or surface ele- 
ment is mapped on a unit line or on a flat unit cell, respectively. 

2. Depending on the proximity between the field point (£) and the element under con- 
sideration, there may be element subdivision and additional mapping for improved 
accuracy. 

3. Gaussian quadrature formulas are employed for the evaluation of the discretized inte- 
gral over each element (or sub-element). These formulas approximate the integral as a 
sum of weighted function values at designated points. The error in the approximation 
is dependent on the order of the (Gauss) points employed in the formula. To mini- 
mize error while at the same time maintaining computational efficiency, optimization 
schemes are used to choose the best number of points for a particular field point and 
element (Watson, 1979). 
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4. When the field point coincides with a node of the element being integrated, the inte- 
gration becomes singular. In this case, the value of the coefficients of the Fij kernel 
corresponding to the singular node cannot be calculated accurately by numerical inte- 
gration. Instead, after the integration of all elements is complete, this value is deter- 
mined so as to satisfy a rigid body displacement of the body (Banerjee and Butterfield, 
1981). 

2.4.3 Assembly of Equations 

After the derivation of the modified boundary integral equations and the analytical 
circumferential integration of the kernel functions, the next critical step in the formulation 
is the assembly of the inserts in the system equations. Here, efficiency is of utmost im- 
portance. The approach to writing an efficient algorithm is to keep the number of system 
equations to a minimum by eliminating all unnecessary unknowns from the system. The 
strategy is to retain in the system only traction variables on the matrix-insert interface. 
This is in contrast to a general multi-region problem where both displacement and trac- 
tions axe retained on an interface. The elimination of the displacement on the interface is 
achieved through a backsubstitution of the insert equations in the system equations which 
are made up exclusively from equations written for the composite matrix (on the outer 
surface and on the surface of the holes). The procedure is described below. 

Equation (2.7) is used to generate a system of equations for nodes on the outer surface 
of the composite matrix and for nodes on the surface of the holes containing the inserts. 
Written in matrix form we have 


On the Matrix Outer Surface: 

G°t° - F°u° + Gt H = 0 

(2.8a) 

On the Matrix Hole Surface: 

G°t° - F°u° + Gt H = Iu H 

(2.86) 


where 

t° and u° are traction and displacement (or flux and temperature) vectors on the outer 
surface of the composite matrix 


13 


t H and u H are traction and displacement (or flux and temperature) vectors on the hole 
I is the identity matrix 

G° and F° matrices contain coefficients from the integration over the outer boundary. 
G matrix contains coefficients integrated about the hole/insert 

Our goal is to eliminate u H from the system. To this end, equation (2.6) is written for every 
node on an insert, collocating slightly outside the boundary of the insert [at a distance of 
(1.25)*(insert radius)] where Cy(£) = 0. 

F I 2 u j = G I 2 t* 

Superscript 12 identifies the equations written at points located slightly outside the bound- 
ary of the inserts. 

Noting u H = u 1 and t H = t 1 we have 

F I 2 u h = -G I2 t H (2.9) 

Post multiplying equation (2.8b) by the F 12 matrix in equation (2.9) yields 

F I2 G°t° - F I 2 F°u° + F I2 Gt H = F I 2 u h (210) 

Equation(2.9) can now be set equal to equation (2.10) and the final form of the system is 
derived. 

On Outer Surface: G°t° - F°u° + Gt H = O 

On Hole: F I2 G°t° - F I2 F°u° + (F I2 G + G I2 )t H = O (2.11) 

At every point on the outer surface, either the traction or the displacement (or the 
temperature or the flux) is specified and on the surface of the hole only the tractions (or 
fluxes) are retained. Therefore, the number of equations in the system are equal to the 
final number of unknowns, and hence, the system can be solved. Thereafter, equation 
(2.8b) is used to determine the displacement on the matrix-insert interface. 
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It should be noted that since the displacement (or temperature) about a particular 
hole is present only in the insert equation corresponding to that hole, backsubstitution can 
be performed one insert at a time in a more efficient manner than backsubstitution of all 
inserts at once. Further, nowhere in the assembly process is a matrix inversion necessary. 
This efficient assembly process was made possible due to the unique formulation of the 
modified boundary integral equations developed earlier in this chapter. 

When the composite matrix is divided into a multi-region model, the above insert 
assembly is performed for each region independently. Thereafter, equilibrium and compat- 
ibility conditions are invoked at common interfaces of the substructured matrix composite. 
After collecting together the known and unknown boundary quantities, the final system 
can be expressed as 

A b x = B b y (2.12) 

where 

x is the vector of unknown variables at outer boundary unknown tractions (or fluxes) 
along the hole/insert interface 

y is the vector of known variables on the outer boundary of the composite matrix, 
and 

A b ,B b are the coefficient matrices 

Standard numerical procedures are used to solve the unknowns in equation (2.12). 
Details are described in the computer development section (Chapter 6). 

Once the unknowns are determined, the resulting displacements (or temperatures) 
on the matrix-insert interfaces can be found using equation(2.8b) rendering a complete 
boundary solution for both the inserts and composite matrix. 

2.5 Interior Quantities 

Once all of the displacements and tractions (or temperature and fluxes) are known on 
the matrix outer surface and on the matrix-insert interface, interior quantities of displace- 


15 


ment, stress and strain (or temperature and flux) can be determined at any point in the 
composite matrix or insert. For displacement (or temperature), either the conventional 
boundary displacement (or temperature) integral equation (2.1) or (2.2) can be employed 
or alternatively the modified equations (2.3) or (2.5) can be used. 

Equations for strains can be derived from the forementioned displacement equations 
and the strain-displacement relations. Thereafter, equations for stress are obtained by 
substituting the resulting strain equations into Hooke’s law. 

The resulting equations, however, are not only invalid on the surface, but also difficult 
to evaluate numerically at points close to it. For points on the surface, the stresses can 
be calculated by constructing a local Cartesian coordinate system with the axes 1 and 
2 directed along the tangential directions and the axis 3 in the direction of the outward 
normal. The stresses referred to these local axes (indicated by overbars) are then given 
by: 


u , Eu (_ . _ \ E _ 

#n = z ^3 + ■; o £ n + e 22 I + z ; — e u 

\-V 1 -U 2 \ / 1 + V 

E 

ai2 = ° 21 = 2(TT^) Cl2 

u - Eu ( _ \ E _ 

922 = ~ iz + \ En + £22 J + T+T C22 

^33 = *3 


( 2 . 13 ) 


cr 32 = CT 23 = ^ 2 


^31 = ^13 = <1 


where E is the Young’s modulus, lij defines the components of the strains in the local axes 
system and U are the traction on the boundary. This method of evaluating the stresses on 
the surface was originally devised by (Rizzo and Shippy, 1968). 
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3. Steady-State Uncoupled Thermoelastic BEM Formulation 

3.1 Introduction 

The boundary element formulation for the steady-state uncoupled thermoelastic anal- 
ysis of perfectly bonded ceramic composite structures is similar to the elastostatic/heat 
conduction BEM formulation of Chapter 2. The primary difference is the inclusion of 
one-way coupling terms which, for the present formulation, requires the solution of the 
heat conduction equation prior to solving the thermoelastic equation set. This, however, 
is favorable since solving two uncoupled subsystems independently is more efficient than 
solving the system as a whole. 


3.2. Boundary Integral Equation Formulation 

Once again, the conventional boundary integral equation is the starting point for the 
uncoupled thermoelastic ceramic formulation. The displacement/temperature for a point 
£ inside the elastic composite matrix is 


Cy(0«<(0 = J [g§(*,0<?(*) - F?(x,Ou°(x) 


dS(x) 


+ E/ [g§(*,o«?(*)-^(*.o«?(*) 

n=l *' S " >■ 


dsr^ix) 


(3.1) 


= 1,2, 3,4 

where 

Gij,Fij are the uncoupled thermoelastic fundamental solutions of the governing differential 
equations of the ceramic matrix of infinite extent [Dargush, 1987] 

Cij are constants determined by the geometry at 4 
Ui,ti are displacements and tractions for i = 1,2,3, and temperatures and fluxes for i = 4 
5,5" are surfaces of the outer boundary of the matrix and the n th hold (left for fiber), 
respectively 

N is the number of individual insert fibers 
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Superscripts O and H identify the quantities on the outer surface of the matrix and the 
quantities on the surface of the hole, respectively. 

The conventional boundary integral equation for displacement /temperature can also 
be written for each of the N insert fibers. For the displacement/temperature at a point £ 
inside the n th insert we can write 

CljiOMO = f gn [g{,(x, 04(*) - F4(x,0«f(z)]dS"(x) (3.2) 

i,j = 1,2, 3, 4 

axe the fundamental solutions of the n th insert 
Cfj axe constants determined by the geometry at £ in insert n 
u{,t{ axe displacement and tractions ( i = 1,2,3) and temperature and fluxes (i = 4) asso- 
ciated with the n th insert 
5" the surface of the n th insert 


We next examine the interface conditions between the composite matrix and the in- 
sert. For a perfect bond the displacement/temperature of the matrix and the displace- 
ment/temperature of the inserts along the interface are equal and the tractions/fluxes are 
equal and opposite. 

uf i (x) = uf(x) (3.3a) 

<f(x) = -tf(x) (3.36) 


Substitution of equations (3.3) into equation (3.2) yields the following modified boundary 
integral equation for insert n. 

cijiOMO = J [~ Gi,'(*,0<?(*) + (*.0“"(*) dsn ( x ) (3- 4 ) 


Finally adding N insert equations (3.4) to equation (3.1) and canceling terms, yields the 
modified boundary integral equation for the composite matrix 


dS(x ) 


CijiOuiiO = J [Gg(x,^)<?(x) - F^(x, 0«?(*) 

n=l' / 5’>L J 


x) 


(3.5) 
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where 


= G» j (x,()-(G I ij (x,0) n ( 3 - 6 ) 

0 = *#(*. 0 + (*«(*. 0 )" ( 3 - 7 ) 

Cy(0 are constants dependent on the geometry for a point ( on the outer boundary and 
Cjj(£) = Sij for a point in the interior of the body. 

3.3 Assembly of Equations 

The approach to writing an efficient algorithm is to keep the number of system equa- 
tions to a minimum by eliminating all unnecessary unknowns from the system. Once 
again the strategy used is to retain in the system only traction/flux variables on the 
matrix-insert interface. This is in contrast to a general multi-region problem where both 
displacement /temperatures, and tractions/fluxes are retained on an interface. The elimi- 
nation of the displacement/temperatures on the interface is achieved through a backsub- 
stitution of the insert equations into the system equations which are made up exclusively 
from equations written for the composite matrix (on the outer surface and on the surface 
of the holes). The procedure is described below. 

Equations (3.4) and (3.5) can be discretized and integrated as described in chapter 2. 
The thermoelastic and heat conduction equations are integrated simultaneously. Therefore, 
using equation (3.4) a system of equations can be generated for the nodal points on the 
discretized composite matrix model. The equations for the nodes on the outer boundary 
of the matrix and for the nodes on the hole surface can be written separately in matrix 
notation as: 


On the Matrix Outer Surface: G°t° - F°u° + Gt H - Fu H = 0 (3.8a) 

On the Matrix Hole Surface: G°t° - F°u° + Gt H - Fu H = Iu H (3.86) 

where 
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t° and u° are traction/flux and displacement/temperature vectors on the outer surface 
of the composite matrix 

t H and u H are traction/flux and displacement /temperature vectors on the hole 
I is the identity matrix 

G° and F° matrices contain coefficients from the integration over the outer boundary. 
G and F matrices contain coefficients integrated about the hole/insert 


Our goal is to eliminate u H from the system. To this end, equation (3.4) is written for every 
node on an insert, collocating slightly outside the boundary of the insert [at a distance of 
(1.25)* (radius of insert)] where C/,(0 = 0. 

G I2 t J - F I2 u j = 0 

Superscript 12 identifies the equations written at points located slightly outside the bound- 
ary of the inserts. 

Noting u H = u 1 , and t H = -t 1 we have 


F I2 u h = -G I2 t H 


(3.9) 


Post multiplying equation (3.8b) by the F 12 matrix in equation (3.9) yields 

F I2 G°t° - F I2 F°u° + F I2 Gt H - F I2 Fu h = F I2 u H (3.10) 

Equation (3.9) can now be set equal to equation (3.10) and the final form of the system is 
derived. 

On Outer Surface: G°t° - F°u° + Gt H — Fu H = O (3.11a) 

On Hole: F I2 G°t° - F I2 F°u° + (F I2 G + G I2 )t H - F I2 Fu H = O (3.116) 

Matrix F of equation (3.11) are coefficients derived from the integration of the Py 
kernel defined in equation (3.7) as 

p . . — pf? pi, 

r *j — r tj t r ij 
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In the elastostatic/heat conduction formulation of Chapter 2 this kernel is zero according 
to equation (2.3c). This allowed the displacement /temperature vector u H to be removed 
from the formulation leaving an equal number of unknowns and equations for which a 
solution can be obtained. In the present formulation, if we once again assume the Poisson 
ratio of the insert to be the same as the composite matrix, all but three of the sixteen 
components in the Fij kernel vanish. F41, F42, and F43 are non-zero when the coefficient of 
thermal expansion of the insert is not the same as the matrix. Examining the Fij kernel 
more closely the following simplification is possible. 


^b( u >O u i (• r ) — 


a Af - a 1 


V M 


■0 

0 

0 

F & 1 



0 

0 

0 

F & 



0 

0 

0 


< 

u i 

.0 

0 

0 

0 J 


lu"J 


(3.12) 


where 


U H U H U H 
) u 2> “3 


,H 




a* 


are the displacement components of the hole 
are the temperature component of the hole 

are components of the F -J kernel on the hole surface of the matrix, and 
are the coefficients of thermal expansion of the insert and the matrix, respec- 
tively. 


From the above equation it is clear the displacement about the hole/insert once again 
vanishes from the formulation, however, the temperature 6 H (6 H = u *£) is present as a 
coupling term, in the displacement equations (components 1,2 and 3 of the kernel) and 
vanishes in the temperature equation (component 4). Equation (3.11) may therefore be 
uncoupled and written as follows. 


Heat conduction system 


Outer Surface: G°t° - F°u° + Gt H = 0 (3.13a) 

Hole Surface: F I2 G°t° - F I2 F°u° + (F I2 G + G I2 )t H = 0 (3.136) 
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Thermal elastic system 


Outer Surface: G°t° - F°u° + Gt H = B0 H (3.13c) 

Hole Surface: F I2 G°t° - F I2 F°u° + (F I2 G + G I2 )t H = B M <? H (3.13d) 


where 


B0 h = Fu H 


and 


B m 0 h = F I2 Fu h 


(3.146) 


In the heat conduction system the number of unknowns is equal to the number of equa- 
tions. The heat conduction system can be solved and the nodal temperatures on the hole 
determined using an uncoupled form of equation (3.8b). With the nodal temperatures on 
the hole known, the number of unknowns in the thermoelastic equation system is reduced 
to exactly the number of equations, and therefore, this system can be solved. With the 
use of equation (3.8b) all remaining nodal boundary variables can be determined on the 
matrix-insert interface rendering a complete solution to the boundary value problem. Dis- 
placement, temperature, traction, flux, stress, or strain measures can now be determined 
at any point on the boundary or in the interior of the body as described in chapter 2. 
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4. Nonlinear Interface Connections Between the Fiber and the Matrix 

4.1 Introduction 

In order to accurately analyze a ceramic composite structure the interaction between 
the fiber and the composite matrix must be properly modeled. Interface phenomena such 
as perfect bonding, progressive debonding (with gap openings), frictional slipping, and 
plastic behavior along the matrix-insert interface must be included. The failure mode of 
an interface is dependent on the state of stress at the interface. Therefore, the general 
mode of failure will be nonlinear and irreversible rendering a path dependent, quasistatic 
analysis. 

In this chapter the steady-state uncoupled thermoelastic boundary element formula- 
tions of chapter 3 is rederived in a form suitable for the inclusion of nonlinear interface 
connections. In addition to the perfectly bonded interface, two new types of interface 
connections are presented in a general form so as to be interchangeable upon input by the 
user. Finally the assembly of the numerically integrated BEM equations is presented and 
an incremental algorithm for their solution is described. The elastostatic and steady-state 
heat conduction formulations are obtained from the uncoupled thermoelastic formulation. 


4.2. Boundary Integral Equation Formulation 

The conventional boundary integral equation for displacement /temperature is the 
starting point for the steady-state uncoupled thermoelastic ceramic formulation with non- 
linear matrix-insert interface connections. The displacement/temperature for a point ( 
inside the elastic composite matrix is 

C^O-iio = J [cg(*. «)«?(*) - fff(x,0«?(*)]<»(*) 

+E/ sn [cg(x,o*r(*)-i^(*.o«i , (*)]-s w (*) (4.i) 

ij = 1,2, 3, 4 

where 
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Gij, Fij axe the fundamental solutions of the governing differential equations of the ceramic 
matrix of infinite extent [Dargush, 1987] 

Cij are constants determined by the geometry at £ 

Ui,u are displacements and tractions 

5, 5" axe the surfaces of the outer boundary of the matrix and the n th hole (left for fiber), 
respectively 

N is the number of individual insert fibers 

Superscripts O and H identify the quantities on the outer surface of the matrix and the 
quantities on the surface of the hole, respectively. 

The conventional boundary integral equation for displacement can also be written for 
each of the N insert fibers. For the displacement/temperature at a point ( inside the n th 
insert we can write 




(4.2) 


i,j = 1,2, 3, 4 

, F/j axe the fundamental solutions based on material properties of the n th insert 
C-j are constants determined by the geometry at £ in insert n 
u{,t( are displacement and tractions associated with the n th insert 
5" the surface of the n th insert 

We next examine the interface conditions between the composite matrix and the insert. 
The difference between the displacement /temperature of the insert and the displacement/ 
temperature of the insert is defined by a variable d%. The traction/flux on the hole is equal 
and opposite to the traction/flux on the insert. Therefore 


uf(z) = «<(*) + <ii(*) 


(4.3a) 

(4.36) 
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Substitution of equations (4.3b) into equation (4.2) yields the following modified boundary 
integral equation for insert n. 




dS"(x) 


(4.4) 


Finally adding N insert equations (4.4) to equation (4.1) and invoking equation (4.3a), 
yields the modified boundary integral equation for the composite matrix 

Cy«K«) + Xj (c$(0)' B «<(0 = j s [<^(*,0<°(*)--F«(*.0“?(*)]^(*) 

+ E / N(x,0^(x)-^(x,0di(x)-^(x,0^(x)]d5"( I ) (4.5) 

n=l Js " J 


where 

e? j (x,() = G? j (x,t)-{G! j {x,0) n 
FT j (x,0 = F${x,0+(F; j (x,0) n 

For a point in the interior of the matrix 

Cy(0 = Sij and (C£( O)" = 0 for all n (4.6a) 

For a point in the interior of the Jfc th insert 

Cij(() = 0 and (Ci«))"={,°. Itl <«»> 

For a point on the outer boundary of the matrix (C§(£) = CV,(0) 

C?(0 is dependent on the geometry at £, and 
(«$<€))" = 0 for all n (4.6c) 

For a point on the k th matrix-insert interface (C^(£) = C^(0) 

c$( 0 +((cfj(< 0) fc = *« 

((Cy(0)" = 0 for ( 4 - 6d ) 


25 


For a point £ on the k th matrix-insert interface the left-hand side of equation (4.5) may be 
rewritten using equation (4.3a) and (4.6d) as 

c(j + (c!j( o) V(o + £ ( c L(0) n ^(0 

n = l 
n*k 

= cjf(e)uf(Q + cg(f)di(t) + (ci(0)V( o = + cgawt) (*v 

4.3 Interface Constitutive Relationships 

4.3.1 Introduction 

The interface constitutive relationships of ceramic composite materials are, in general, 
nonlinear. This requires that their BEM solution be found using an incremental (qua- 
sistatic) algorithm. Therefore, the linear BEM equations of Section 4.2 are interpreted as 
incremental relations and the interface constitutive relationships are defined in terms of 
the incremental displacements/temperatures Uj and the incremental tractions/fluxes £, as 

U = kpj (4.8a) 

i,j = 1 to 3 for elastostatics 
i,j = 1 for heat conduction 
i,j = 1 to 4 for uncoupled thermoelasticity 

where k*? is the nonlinear constitutive matrix dependent on the current state of stress, 

dj = uf - uj, (4.86) 

U ~ t\ = -ii (4.8c) 

The elastostatic constitutive relationships that are derived in this section are expressed 
using a local coordinate system where the first component corresponds to the direction nor- 
mal to the interface and the second and third components correspond to arbitrary tangen- 
tial directions on the matrix-insert interface. The local constitutive matrix is transformed 
to the global coordinate system when incorporated in the BEM equations as follows 
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where a mj is the directional cosine transformation tensor. 


4.3.2 Linear Spring Interface Connection 

The general form of the linear spring interface connection may be expressed as (kfj 
replaces kf?) 

U = kM (4-9) 

In the current elastostatic implementation the special case k{^ = k^k^ = fcfa = k t , and k £ = 0 
for i ? j, is assumed (see Figure 4.1). k n is the spring constant in the normal direction 
and kt is the spring constant in the tangential direction. The general form does not pose 
any more difficulty in implementation then the special case, but the additional generality 
is usually not required in practice. The limiting form as k n —* oo and k t -* oo approaches 
the case of a perfect bond, and k n — k t = 0 corresponds to a completely debonded interface. 
jfc n — + oo and k t = 0 corresponds to a sliding interface with a perfect connection in the normal 
direction. Results using the present formulation with values k n = k t = (£' matnx + E"" 8 * 1 *) * 10 6 
compared very closely to results obtained for the same problem assuming a perfect bond. 
Similar comparisons were found to hold true for problems with k n — k t = 0 verses £ insert = o 
(hole solution). 

In the heat conduction implementation fcn = —k where k is the thermal conductivity, 
it = 0 does not permit heat flow across the interface which is analogous to am insulated hole 
and Jb -» oo is analogous to a perfect bond between the insert and the matrix. 

The kij for uncoupled thermoelastic analysis is derived by combining the elastostatic 
and heat conduction kij to form a completely uncoupled k^ matrix k n = k n ,k ^2 = *33 = 
Jb t , Jfc 44 = -Jfc, and k^ = 0 for i ± j. 

4.3.3 Spring-Friction Nonlinear Interface (Coulomb Friction) 

An interface model shown in Figure (4.2) exhibits a linear spring behavior normal to 
the interface when the normal tractions are in compression. Linear spring resistance is 
also observed in the tangential direction when the principal tangential traction is below 
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the slip limit defined by the coulomb friction criteria. When the slip limit is reached the 
traction that the interface can sustain in the tangential direction has reached a maximum 
and any additional tangential traction will result in an irrecoverable shift of the relative dis- 
placement between the insert and the matrix, and a redistribution of stress (and therefore 
interface tractions) is required to bring the structure back in equilibrium. Furthermore, 
the model does not support tensile tractions normal to the interface. Instead the tractions 
at this point are set to zero through the constitutive relation kij = 0 (i,j — 1 to 3) and a 
gap between the insert and the matrix will form. Once again a redistribution in stress is 
required to bring the structure back in equilibrium. 

The nonlinear interface constitutive relationship for a point on the interface exhibit- 
ing the spring- coulomb friction phenomena can be derived in a manner analogous to the 
incremental theory of plasticity [Selvadurai, 1988]. This requires a description and use of 

(1) A flow rule relating the irrecoverable nonlinear part of the displacement difference rate 
to the traction state at the interface. 

(2) A consistency relation requiring the new traction state to lie on the newly formed yield 
surface (defined by the hardening rule and yield function). 

(3) A yield function defining the limits of elastic behavior, and 

(4) A hardening rule defining the subsequent yield surfaces. 

The variable for the incremental displacement difference across the interface is assumed 
to be composed of an elastic part tf? and a plastic part d? 

d^dt + d? (4.10) 

where the elastic component is related to the interface traction by the linear constitutive 
relation defined in Section 4.3.2. 

i . - ]c^ d e 
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or 


Next a flow rule is defined as 


U = %&-%) 


i-ig 


( 4 . 11 ) 


J atj 

where A is an unknown flow factor dependent on the current state of traction and Q is a 
nonlinear potential. Substituting the flow rule into equation (4.11) yields 


t - — k e d - — k e 

i, _ Ky-a j k 'j Qij 


( 4 . 12 ) 


The consistency relation is defined as 


dF. n 

— — ti = 0 
dti 


( 4 . 13 ) 


Substitution of equation (4.12) into the consistency relationship and rearranging yields 
a relationship for A. 

( 4 . 14 ) 


x ~ Gdti J 


where G = m^k' mn §2 


Finally substitution of equation (4.14) into equation (4.12) yields 


ti — kijdj 


( 4 . 15 ) 


where 


and 


jfcf? - jt«, - jb?. 
K iy ~ K ij 


k ? -L k ' k 

K ij- G * p* dtp dt q 


( 4 . 16 ) 


( 4 . 17 ) 


The nonlinear interface constitutive relationship is complete once F and Q are defined. 
Assuming the coulomb friction criteria, the yield function F is given by 


F — (<2 + tl) 1 ^ 2 + tdn — 0 


( 4 . 18 ) 
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where < n is the total traction normal to the interface, * 2 and <3 are the total tangential 
tractions and /i is the coefficient of friction between the insert and the matrix. 

This relation holds only for t n < 0. If F < 0, the interface is assumed elastic and a 
spring connection is used. If F = 0 the nonlinear interface relation (4.15) is used. F < 0 
is meaningless in theory, however, since in practice we assume a finite size load step the 
function may become slightly positive in which case it is treated as F = 0. When t n > 0, 
kij = 0 (i,j = 1 to 3) is assumed. 

This relationship is assumed constant for the entire analysis and therefore to complete 
the analogy with classical plasticity, we can say we have zero hardening. The nonlinear 
potential Q is assumed to be the magnitude of the principal tangential traction. 

Q = (tl + t 2 3 y' 2 (4.19) 


If we assume the special elastostatic case described in Section 4.3.2 for A- where An — 
k n , A 22 = A 33 = A t , and A?- = 0 for i ^ j we can write for kF 


jfcP - 
K ij 


0 0 0 
fik n t2 k it 2^2 A t <2<3 

A n 1 3 A £^ 2 1 3 At* 3*3 


(4.20) 


where <2 and <3 axe the normalized components of total tangential tractions, i.e. 


<2 = 


1 2 


and 


<3 = 


(*!+*!) 1/2 


In the present heat conduction and uncoupled thermoelastic implementation, only a lin- 
ear (uncoupled) relation is allowed between the increment of flux and the temperature 
difference rate variable across the interface. This of course will be updated in the newly 
proposed work. 


4.4 Assembly of Equations for General Matrix-Insert Interface Connections 
In the assembly of the perfectly bonded inserts the number of equations in the final 
system was reduced by eliminating all unnecessary unknowns from the problem. In the 
previous case all unknown displacement / temperature nodal variables on the matrix-insert 
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interface were eliminated leaving only the unknown tractions/fluxes. In the present case the 
tractions/fluxes will be eliminated retaining only the matrix- insert displacement difference 
variables along the interface. 

The steady-state uncoupled thermoelastic equations (4.4) and (4.5) can be discretized 
and integrated as described in Chapter 2. Therefore, using equations (4.4), (4.6c), and 
(4.7) a system of equations can be generated for the nodes on the discretized model of the 
composite matrix. The equations for the nodes on the outer boundary of the matrix and 
for the nodes on the hole surface can be written separately in matrix notation as: 

On the Matrix Outer Surface: G°t° - F°u° 4- Gt H - F H d - Fu 1 = 0 (4.21a) 

On the Matrix Hole Surface: G°t° — F°u° + Gt H - F H d - Fu 1 = Iu 1 (4.216) 

where 

t° and u° are traction/flux and displacement/temperature vectors on the outer surface 
of the composite matrix 

t H is the traction/flux vector on the hole 
u 1 is the displacement/temperature vector on the insert 
d vector of displacement /temperature difference between the hole and insert 
I is the identity matrix 

G° and F° matrices contain coefficients from the integration over the outer boundary 
G, F and F H matrices contain coefficients from the integration over the hole/insert 

Vectors u° and d are arranged in the same order as the equations are written so as to 
produce a strong main diagonal in matrices F° and F H . 

Equation (4.6c) is used in writing equation (4.21a). The cQ(()u ° (£) term is lumped with 
its respective (diagonal block) coefficient in matrix F°. The sum of these two terms are 
accurately calculated indirectly using the rigid body technique [Banerjee and Butterfield, 
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1981]. Equation (4.7) is used in writing equation (4.21b). The C{j(£)di{£) term is lumped 
with its respective (diagonal block) coefficient in matrix F H (also calculated by the rigid 
body technique). Similarly, the term JyuftO could be lumped with matrix F, but instead 
it is put on the right-hand side (forming an identity matrix I for use in subsequent matrix 
algebraic operations. Next, equation (4.5) is written for every shape function node on an 
insert, collocating slightly outside the boundary of the insert [at a distance of (1.25)*(radius 
of insert)] where Cf (£) = 0. 

G^t 1 - F I 2 u ! = 0 (4.22) 

Superscript 12 identifies the equations written at points located slightly outside the bound- 
ary of the inserts. Equation (4.22) can be rewritten as 

F I 2 u i = -G I2 t H (4.23) 

Post multiplying equation (4.21b) by matrix F 12 of equation (4.23) yields 

F I2 G°t° - F I2 F°u° + F I2 Gt H - F I2 F H d - F I2 Fu* = F I2 u* (4-24) 

For efficiency, this operation can be carried out one insert at a time. Equation (4.23) can 
now be set equal to equation (4.10). This yields, 

F I2 G°t° - F 12 F°u° + G M t H - F M d - F^Fu 1 = 0 (4 25) 


where 

g m _ F I2 G+ G 12 , and 

F M _ p I 2 p H 

The nonlinear interface relation t H = -K ep d is substituted in equation (4.21a) and (4.25) 
and noting equations (3.12) and (3.14) the final form of the system can be written as 

On Outer Surface: G°t° - F°u° - (GK ep + F H )d = B0 1 (4.26a) 

On Hole Surface: F I2 G°t° - F I2 F°u° - (G M K ep + F M )d = B M ^ (4.266) 
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For each degree of freedom on every node on the outer surface, either the displace- 
ment/temperature or traction/flux is specified and on the matrix-insert interface the dis- 
placement/temperature difference variable and the insert temperature variable are present. 
Therefore the number of unknowns in the system is greater than the number of equations 
by exactly the number of nodes on the interface (exceeded due to the insert temperature 
variables). However, as was the case in uncoupled thermoelasticity with a perfectly bonded 
matrix-insert interface, the B and B M matrices influence only the thermoelastic solution, 
i.e., the coefficients corresponding to the heat conduction equations are zero. Therefore, in 
the heat conduction equations these terms do not appear. The heat conduction equation 
can be uncoupled from the system in a manner similar to equation (3.13) rendering an 
equal number of heat conduction equations and unknowns. The heat conduction system is 
solved and the insert nodal temperatures are determined using equation (4.21b). With the 
insert nodal temperatures known, the number of unknown in the thermoelastic equation 
system is reduced to exactly the number of equations, and therefore the system can be 
solved. With the use of equations (4.8) and (4.21b) all remaining nodal boundary vari- 
ables can be determined on the matrix-insert interface rendering a complete solution to 
the boundary value problem. Displacement, temperature, traction, flux, stress, or strain 
measures can now be determined at any point on the boundary or in the interior of the 
body as described in Chapter 2. 

The elastostatic and steady-state heat conduction formulations are derived from equa- 
tion (4.26). Noting the F^ kernel vanishes in elastostatic and heat conduction analysis, 
equation (4.21) is used setting matrix F = 0. The final system is similar to equation (4.26) 
with B = B m = 0. The number of unknowns are equal to the number of equations in a well 
posed problem and the system can be solved in one step. 

4.5 Numerical Algorithm for Nonlinear Interface Connections 

The nonlinear constitutive relationships are incorporated in the linear BEM equations 
and solved using the incremental (quasistatic) algorithm described below. 
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Nonlinear Solution Algorithm 

(a) The boundary (and body force loading if present) is divided into a number of small 
sub-increments. Ten sub-increments have been found to be sufficient. 

Initial Load Step (Iterative Start-up Loop): 

(b) Assemble the BEM system equation (4.26) as described in Section 4.4 using the in- 
terface relations calculated in step (d). In the initial assembly step, linear spring 
connections (Section 4.3.2) are assumed at each interface node. 

(c) Apply the sub-increment of load to the system equations (4.26). 

(d) Determine the traction state at each interface node and derive the appropriate interface 
constitutive relation for that node. 

- If the normal traction a the interface is positive, assume a gap opening with zero 

traction, i.e., = 0 

- If F < 0 assume a spring connection (Section 4.3.2) 

- If F < 0 assume a nonlinear connection (Section 4.3.3) and evaluate the constitutive 
relationship using the current traction state (from current loading). 

(e) Changing the interface connections from a linear spring to nonlinear relationship in the 
first sub-load step may cause a major redistribution in the tractions on the interface. 
Therefore, the constitutive relations calculated in step (c) will change. Repeat step 

(b) , (c) and (d) until a coverage solution is achieved for the first load step. 

Subsequent Load Steps: 

(f) Steps (b), (c), and (d) are repeated once for each additional load step. In general, 
iteration within a subsequent sub-load step is not required since the constitutive re- 
lationships at interface nodes change in a gradual manner once the start-up loop is 
complete. The solution of any sub-increment is the accumulation of all previous sub- 
load steps. It is this current solution that is used in step (d) to evaluate the current 
interface constitutive relationships used in the next sub-incremental load step (step 

(c) ). 
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5. TRANSIENT HEAT CONDUCTION AND TRANSIENT UNCOUPLED 
THERMOELASTIC BEM FORMULATIONS 

5.1 Introduction 

The analysis of a body in a transient state is inherently more complicated than a 
steady-state analysis. The boundary element formulation for transient analysis contains 
convoluted integrals which must be integrated in time as well as in space. Therefore, a 
transient BEM analysis is expected to be more complicated and more expensive than its 
steady-state counterpart. To compound these difficulties, the cost saving techniques devel- 
oped for the steady-state ceramic composite analyses cannot be employed for the transient 
case. First, the efficient modified boundary integral formulations and the resultant as- 
sembly schemes of Chapters 2-4 do not lend themselves to transient insert formulations 
for composites with general material properties. Secondly, due to the complexity of the 
transient kernel functions, the integration about the circumference of the insert must now 
be carried out using numerical integration, thus adding to the cost of the analysis. 

Nevertheless, an efficient numerical integration scheme as well as an efficient assembly 
algorithm have been derived for the transient heat conduction and thermoelastic BEM 
analyses. In the assembly algorithm the equations of the inserts are backsubstituted into 
the equation system of the composite matrix, hence, reducing the size of the overall system. 
This is carried out prior to decomposition at a point when the equations for each insert 
can be handled individually in the most efficient manner. This assembly, as well as the 
decomposition, is required only for the initial time step. Later time steps require only a 
new calculation of the system’s right-hand side. 

5.2 Transient Boundary Integral Equation Formulation 

The transient, uncoupled thermoelastic boundary integral equation (Daxgush and Baner- 
jee, 1991) for the displacement (and temperature) at a point £ in a composite matrix is 
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( 5 . 1 ) 


C 0 a ( , 0 = / be * #(X, 0 - /*, * 1#(X, <)] dS(X) 

JS 
M 

+ E / to** * *"(*. 0 - fpc * u» (X , <)] d<r(x) 

m = 1 • /5m 

where 

a,/? indices varying from 1 to 4 for uncoupled thermoelasticity and equaled to 4 for 
heat conduction 

Cp a constants determined by the geometry at ( 

■up,tp generalized displacement and traction 
Up = [ui, u 2 , u 3 , 6} T 
tp=[t 1 , < 2 , < 3 , q] T 
0,q temperature, heat flux 

gap, fap generalized displacement and traction kernels 
S, S'” are the surfaces of the outer boundary of the matrix and the mth hole (left for 
fiber) 

M is the number of individual fibers 
and, for example 

90a * tp = [ 9pa{X, t-,Z, T)tp(X , r)dr 

Jo 

denotes a Riemann convolution integral. 

Superscripts O and H are used to highlight the quantities associated with the outer 
surface of the matrix and the quantities associated with the surface of the hole, respectively. 
The boundary integral equation for mth insert can be written as 

= [ [900 * 4(X, t) - f ! p a * u’p(X, <)] dS™(X) ( 5 . 2 ) 

J S m 

where 

Qpa’fpa the fundamental solutions for a body with material properties corresponding 
to the mth insert. 
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C/jo axe constants determined by the geometry at £ in insert m 
Up, t r p axe the generalized displacement and tractions associated with the mth insert 

S™ the surface of the mth insert 

In principle, at each instant of time progressing from time zero, these equations can be 
written at every point on the respective boundary. The collection of the resulting equations 
could then be solved simultaneously, producing exact values for all the unknown boundary 
quantities. In reality, of course, discretization is needed to limit this process to a finite 
number of equations and unknowns. Techniques useful for the discretization of (5.1) and 
(5.2) are the subject of the following section. 

5.3 Numerical Implementation 

The boundary integral equations (5.1) and (5.2) developed in the last section, are exact 
statements. No approximations have been introduced other than those used to formulate 
the boundary value problem. However, in order to apply (5.1) and (5.2) for the solution 
of practical engineering problems, approximations are required in both time and space. In 
this section, an overview of a general-purpose, state-of-the-art numerical implementation 
is presented for transient analysis. 

5.3.1 Temporal Discretization 

Consider, first, the time integrals represented in (5.1) and (5.2) as convolutions. Clearly, 
without any loss of precision, the time interval from zero to t can be divided into N equal 
increments of duration at At. 

By assuming that the primary field variables, tp and up, are constant within each At 
time increment, these quantities can be brought outside of the time integral. That is, 

N *nAt 

gpa*tfiiX,t) = J2q{X) / 9 p a (X -£,t - r)dr (5.3a) 

n=l J 

N rn&t 

fpa * MX, t) = T u n p{X) / fMX -t,t-T)dr (5.36) 

nTl J(n- l)At 
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where the superscript on the generalized tractions and displacements represents the time 
increment number. Notice, also, that within an increment, these primary field variables 
are now only functions of position. Since the integrands remaining in (5.3) are known in 
explicit form from the fundamental solutions, the required temporal integration can be 
performed analytically, and written as 


rr\At 

G" +1 -"(X-0= / 9fia(X-t,t-T)d7 

J (n— l)At 


(5.4a) 


(5.46) 


rni it 

F p a +1 ~ n (x - 0 = / /*»(* - * . * - T ) dr - 

These kernel functions, G$JX - 0 and Fg a (X-Q, are detailed in Appendix B. Combining 
(5.3) and (5.4) with (5.1) and (5.2) produces 

Cp a ( 0^(0 = / [Gg 1 ^{x-t)tftX)-F™- n {x-t)uftxj\ dS(X) 

n= 1 ^ Js 

+ E / te +1_ "(^ - OW - F p a +1 ~ n ( X - 0«3(*)] dSm (X) } (5.5 a) 

m=l ds m > 


and for the mth insert 


C*,(0«?(0 = E { / m [ G "a 1_ "(* - * WO - F p a +1 ~ n ( x - 0«2(*)] d!T{X) } (5.56) 

n=l ^ J S m 

which are the boundary integral statements after the application of the temporal discretiza- 
tion. 


5.3.2 Spatial Discretization 

With the use of generalized primary variables and the incorporation of a piecewise 
constant time stepping algorithm, the boundary integral equation (5.5) begins to show a 
strong resemblance to the equations of the steady-state analyses, particularly for the initial 
time step (i.e., N = l). In this subsection, those similarities will be exploited to develop 
the spatial discretization for the uncoupled quasistatic problem. This approximate spatial 
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representation will, subsequently, permit numerical evaluation of the surface integrals ap- 
pearing in (5.5). The techniques described here, actually, originated in the finite element 
literature, but were later applied to boundary elements by Lachat and Watson (1976). 

The process begins by subdividing the entire surface of the body (including the inserts) 
into individual elements of relatively simple shapes. The geometry of each element is, then, 
completely defined by the coordinates of the nodal points and associated interpolation 
functions, as described in Chapter 2. 

In the present work, the geometry is exclusively defined by quadratic shape functions. 
On the other hand, the variation of the primary quantities can be described, within an ele- 
ment, by either quadratic or linear shape functions. (The introduction of linear variations 
proves computationally advantageous in some instances.) 

Once this spatial discretization has been accomplished and the body has been subdi- 
vided into Q boundary elements and P insert elements, the boundary integral equation can 
be rewritten for the matrix as 

E [ / G^- n (X(0-0UC)dSHX(0)] t%, 

,=1 Vs' J 

- E [f S ' F 0a +1 ~ n WO - OMCMSWO)] ^ 

+ E [ / G "a +1 “"(*(0 - 0NM)M y (6)d<?(X(C), 0)1 

p= 1 Usr J 

- E \l W'-ixiO -OWWW),!)] (5.6a) 

and 

E [ / G £ +1_r W) -0^(0^W^(^(0.0l y 

- E [J SP -0NU0M-r(e)dSP(X(0,6^ (5.66) 

where 


CpdOupiO 


N r 

n= 1 ^ 




N r 

n=l '• 
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Q is the number of boundary elements on the outer surface of the composite matrix 
in the region, 

P is the number of insert elements in the region, 
represents a two-dimensional shape function, 

N 7 represents a one- dimensional shape function, 

M y is the circular shape function defined in Chapter 2, 

t,u are nodal values of generalized traction and displacement, respectively. 

In the above equation, the nodal quantities and were brought outside the surface 
integrals. This positioning of the nodal primary variables outside the integrals is, of course, 
a key step since now the integrands contain only known functions. However, before dis- 
cussing the techniques used to numerically evaluate these integrals, a brief discussion of 
the singularities present in the kernels Gp a and Fp a is in order. 

The fundamental solutions to the uncoupled quasistatic problem contain singularities 
when the load point and field point coincide, that is, is when r = 0. The same is true of 
and Fp a , since these kernels are derived directly from the fundamental solutions. Series 
expansions of terms present in the evolution functions can be used to deduce the level of 
singularities existing in the kernels. 

A number of observations concerning the results of these expansions should be men- 
tioned. First, as would be expected has a stronger level of singularity than does the 
corresponding G l a p, since an additional derivative is involved in obtaining F^p from G^p. 
Second, the coupling terms do not have as high degree of singularity as do the correspond- 
ing non-coupling terms. Third, all of the kernel functions for the first time step could 
actually be rewritten as a sum of steady-state and transient components. That is, 

G 1 S3/~v . tr/~r 1 

a /3 ~ U a /3 + ap 

771I as T? 1 tr zpl 

tc*P + t a (3 

Then, the singularity is completely contained in the steady-state portion. Furthermore, 
the singularity in G]j and F^ is precisely equal to that for elastostatics, while G\ e and 
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singularities are identical to those for potential flow. This observation is critical in the 
numerical integration of the F a p kernel to be discussed in the next subsection. However, 
from a physical standpoint, this means that, at any time t, the closer one moves toward the 
load point, the closer the quasistatic response field corresponds with a steady-state field. 
Eventually, when the sampling and load points coincide, the quasistatic and steady-state 
responses are indistinguishable. As a final item, after careful examination of Appendix 
B, it is evident that the steady-state components in the kernels and with n > 1, 
vanish. In that case, all that remains is a transient portion that contains no singularities. 
Thus, all singularities reside in the 33 G a p and 33 F a p components of G l a0 and respectively. 

5.3.3 Numerical Integration 

Having clarified the potential singularities present in the coupled kernels, it is now 
possible to consider the evaluation of the integrals in equation (5.6). 

To assist in this endeavor, the following three distinct categories can be identified. 

(1) The point £ does not lie on the element m. 

(2) The point £ lies on the element m, but only non-singular or weakly singular integrals 
are involved. 

(3) The point lies on the element m, and the integral is strongly singular. 

In practical problems involving many elements, it is evident that most of the integration 
occurring in equation (5.6) will be of the category (1) variety. In this case, the integrand 
is always non-singular, and standard Gaussian quadrature formulas can be employed. So- 
phisticated error control routines are needed, however, to minimize the computational 
effort for a certain level of accuracy. This non-singular integration is the most expensive 
part of a boundary element analysis, and, consequently, must be optimized to achieve an 
efficient solution. In the present implementation, error estimates, based upon the work of 
Stroud and Secrest (1966), are employed to automatically select the proper order of the 
quadrature rule. Additionally, to improve accuracy in a cost-effective manner, a graded 
subdivision of the element is incorporated, especially when £ is nearby. 
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The integration over the surface of the hole and insert elements must be carried out 
using numerical integration. The complexity of the transient kernel function prohibits an 
analytic integration about the circumference of the insert in a manner similar to the steady- 
state case. Therefore, numerical integration must be performed in both the longitudinal 
and the circumferential directions. An efficient integration scheme is adopted in which one- 
dimensional Gaussian integration formulas are applied independently in the two directions. 
This allows the subsegmentation and mapping, as well as the order of the Gaussian formulas 
to vary independently in the two directions so that the accuracy and the cost of the 
integration can be optimized. 

Turning next to category (2), one finds that again Gaussian quadrature is applicable, 
however, a somewhat modified scheme must be utilized to evaluate the weakly singular 
integrals. This is accomplished through element subsegmentation about the singular point 
so that the product of shape function, Jacobian and kernel remains well behaved. 

Unfortunately, the remaining strongly singular integrals of category (3) exist only in 
the Cauchy principal value sense and cannot, in general, be evaluated numerically, with 
sufficient precision. It should be noted that this apparent stumbling block is limited to the 
strongly singular portions, ”Fij and 44 Fee, of the kernel. The remainder of including 
tr Fij and tr Fg e , can be computed using the procedures outlined for category (2). However, 
as will be discussed in the next subsection, even category (3) “Fij and 44 Fee kernels can be 
accurately determined by employing an indirect ‘rigid body’ method originally developed 
by Cruse (1974). 

5.3.4 Assembly of Equations 

The complete discretization of the boundary integral equation, in both time and space, 
has been described, along with the techniques required for numerical integration of the 
kernels. Now, a system of algebraic equations can be developed to permit the approximate 
solution of the original quasistatic problem. This is accomplished by systematically writing 
(5.6) at each node on the outer surface, on the surface of the hole, and on the surface of 
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the insert. The ensuing nodal collocation process, then produces a global set of equations 
given here in matrix form, for time N. 


For the nodes on matrix: 

[<&{#} - [*£]{*£} + [G^H#} - [Fkl{u n „} = (5.7a) 

For the nodes on the inserts: 

[G)}{tY}-[Fj}{u?} = {D?} (5.7 6) 

where {D N } and {D^} are the values of the convolution from the previous N - 1 time steps 
defined as 

~{D n } = £ {(GQ +1-n ]{<2)} - [F£ +1 -"]{ U S} + [G^ +1 -"]m - [F£ +1 -"]{u£}} (5.8a) 

n — 1 

and 

-{D?} = X? {[Gf +1 -"]{<?} - [Ff +1 -"]{a?}} (5.86) 

n= 1 

In these equations, subscripts O, H and I are used to denote the quantities associated with 
the outer surface of the matrix, the surface of the hole containing the insert, and surface 
of the insert, respectively. Furthermore, the C a p(() term of each equation is lumped with 
the respective coefficient in the F 1 matrix. 

To calculate the coefficients of the singular points by the ‘rigid body’ technique, con- 
sider now, the first step. Thus, for N = 1, equation(5.7) becomes 

[G'oUtb) - [F&KuM + [GkHtM - = {0} (5.9a) 

[G}]{4}-[F / 1 ]{u}} = {0). (5.96) 

At this point, the coefficients of the F matrices corresponding to the singular node of the 
equations has not been completely determined due to the strongly singular nature of the 
kernel function. To determine the values of these coefficients, we first decompose the F^ 
kernel into transient and steady-state parts. Following Cruse (1974) and, later Banerjee 
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et al (1986) in elastodynamics, the steady-state part of the coefficients can be calculated 
indirectly by imposing a uniform ‘rigid body’ generalized displacement field on the same 
body under steady-state conditions. The steady-state part of the singular coefficients 
are simply the summation of the non-singular, steady-state coefficients. The remaining 
transient portion of the coefficients are non-singular, and hence can be evaluated to any 
desired precision. 

For an insert that is perfectly bonded to the matrix, the displacement of the surface 
on the hole containing the insert is equal to the displacement of the surface of the insert, 
and the surface tractions are equal in magnitude but opposite in direction. 


{u n H } = {u?} (5.10a) 

{<&} = -{<?} ( 510fc ) 

Substitution of equation (5.10) into equations (5.7b) and (5.8b) and rearranging yields 

{<£} = -[#][*/]{“£} - MP"} ( 5 - lla ) 

{D?} = £ + [F," +1 -"]{u&}} (5.116) 


W ) = [ G }]- 1 

Note, since the boundary integral equation for each insert is independent of the equations 
for the other inserts, the inversion of [G)] can be reduced to a number of smaller inversions, 
one corresponding to each insert. 

Equation (5.11a) can now be back substituted in equation (5.7a) to yield 

[GJ,]{C } - [FoH«o} - [F M ]{u£} = {D»} + [G},][J7KZ?n (5.12) 

where 

[F m ] = [G},][F][F/] + [Fh] 
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In a well-posed problem, at time At, the set of global generalized nodal displacements 
and tractions will contain exactly (4 x P) unknown components (and P unknowns for heat 
conduction analysis). 

Then, as the final stage in the assembly process, equation (5.12) can be arranged to 
form 

= [B l ]{y N ] + {D N } + [GttmDTY ( 5 - 13 ) 

in which 

{x N } unknown components of {u N } and {t^} 

{y^} known components of {u N } and {t^} 

[A 1 ,]!# 1 ] associated matrices 

Note, the entire matrix [F M ] becomes part of [A 1 ] since all the (generalized) displacements 
on the interface are unknown quantities. 

5.3.5 Solution 

A solution of equation (5.13) may be achieved for any time using a time marching 
algorithm. For the initial time N = 1 , {D 1 } = {D l H } — {0} and equation (5.13) can be 
rewritten as 

[A 1 ]** 1 } = [B'Hy 1 } (5.14) 

To obtain a solution for the unknown nodal quantities of this equation, a decomposition 
of matrix [A 1 ] is required. In general, [A 1 ] is a densely populated, unsymmetric matrix. 
The out-of-core solver, utilized here, was developed originally for elastostatics from the 
LINPACK software package (Dongarra et al, 1979) and operates on a submatrix level. 
Within each submatrix, Gaussian elimination with single pivoting reduces the block to 
upper triangular form. The final decomposed form of [A 1 ] is stored in a direct-access file 
for reuse in subsequent time steps. Backsubstitution then completes the determination of 
{x 1 }. Additional information on this solver is available in Banerjee et al (1985). 

After returning from the solver routines, the entire nodal response vectors, {u 1 } and 
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{t 1 }, at time At are known. For solutions at later times, a simple marching algorithm is 
employed. Assuming that the same set of nodal components are unknown as in (5.14) for 
the first time step, equation (5.13) is reformulated as 

[A 1 ]^ 2 } = [B 1 ]*!/ 2 } + {Z? 2 } + (5.15) 

Since the right-hand side contains only known quantities, (5.15) can be solved for {x 2 }. 
However, the decomposed form of [A 1 ] already exists on a direct-access file, so only the 
relatively inexpensive backsubstitution phase is required for the solution. 

The generalization of (5.15) to any time step N is simply equation (5.13) 

[A 1 ]!*"} = [B 1 ]{y JV } + {D N } + [G),][tf]{Z?n (5.16) 

in which the vectors { } and { contain summations which represent the effect of past 
events. By systematically storing all of the matrices and nodal response vectors computed 
during the marching process, surprisingly little computing time is required at each new 
time step. In fact, for any time step beyond the first, the only major computational task 
is the integration needed for form [G w ] and [F^]. Even this process is somewhat simplified, 
since now the kernels are non-singular. Also, as time marches on, the effect of events that 
occurred during the first time step diminishes. Consequently, the terms containing [G N ] 
and [F N ] will eventually become insignificant compared to those associated with recent 
events. Once that point is reached, further integration is unnecessary, and a significant 
reduction in the computing effort per time step can be achieved. 

It should be emphasized that the entire boundary element method developed, in this 
section, has involved surface quantities exclusively. A complete solution to the well-posed 
linear uncoupled quasistatic problem with composite inserts can be obtained in terms 
of the nodal response vectors, without the need for any volume discretization. In many 
practical situations, however, additional information, such as, the temperature at interior 
locations or the stress at points on the boundary, is required. The next section discusses 
the calculations of these quantities. 
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5.4 Interior Quantities 

Once equation (5.16) is solved, at any time step, the complete set of nodal displace- 
ments and tractions are known. Subsequently, the response at points within the body can 
be calculated in a straightforward manner. For any point £ in the interior, the generalized 
displacement can be determined from (5.6a) or (5.6b) with Cp a = 6p„. Now, all the nodal 
variables on the right-hand side are known, and, as long as £ is not on the boundary nor 
on the interface between the insert and the matrix. The kernel functions in (5.6) remain 
non-singular. However, when £ is on the boundary or on the interface, the strong singular- 
ity in ** Fp a prohibits accurate evaluation of the generalized displacement via (5.6), and an 
alternate approach is required. The apparent dilemma is easily resolved by recalling that 
the variation of surface quantities is completely defined by the elemental shape functions. 
Thus, for points on the outer boundary of the matrix, the desired relationship is simply 

«?(© = ( 517a ) 

Where L u are the shape functions for the appropriate element and ( are the intrinsic coor- 
dinates corresponding to ( within that element. Obviously, from (5.17), neither integration 
nor the explicit contribution of past events are needed to evaluate generalized boundary 
displacements. 

In many problems, additional quantities, such as heat flux and stress, are also impor- 
tant. The boundary integral equation for heat flux, can be written 

?r( o = E { E f / £ j&i 1_n wo - oMcyswo)] t% 

n=l ^ q=\ '■•'S' J 

- E [J s , - OMOrfsWO)] «&, 

+£[/, E peV "(X(C)-ON u ,(C)M Y W dsP ( j,f (C))| */L-y 

- f; [J sr D% e f- n (X(0 - 0M0M y (8)dSP(X(0)\ (5.18) 
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where 


wo-o=~* 


dG n 0e (X( 0-0 

% 


(5.19a) 

(5.196) 


A similar equation can be written for the heat flux inside an insert. This is valid for interior 
points, whereas, when £ is on the boundary, the shape functions can again be used. In this 
latter case, 


(0 

3Iw( 0 nN _ ^® Xi n N(c\ 


(5.20a) 

(5.206) 


8( ' w ' k d( 

which can be solved for boundary flux. Meanwhile, interior stresses can be evaluated from 


rff(0 = E { E [ f £ N+1_r W) - OMO^WO) 

n= 1 ^ q= 1 L*' 5 ’ 

- E [ / Dft 1 - - n (x(0 - 0^(C)^(A(C))1 

9=1 ‘■•'S* 

+ E [ / ^-(.(O - t)N»{C)W)d£?{x(Q) tpUTf 

P =1 ^ JSP 

- E [/ SP - o^(c)m 7 (^sp(x(o)] 


(5.21) 


in which 


2/ii/ c 

L - 2i/ <J ‘ du 



(?* 

% J 

| ~~ P$ij G'pe 

(5.22a) 

+ H 

(dF» 
V Kj 

dF eA 
db 1 

— P&ijFpe- 

(5.226) 


A similar equation could be derived if stress measurements are required inside an insert. 
Equation (5.22) is, of course, developed from (5.6). Since strong kernel singularities appear 
when (5.22) is written for outer boundary points, an alternate procedure is needed to 
determine surface stress. This alternate scheme exploits the interrelationships between 
generalized displacement, traction, and stress and is the straightforward extension of the 
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technique typically used in elastostatic implementation (Dargush and Banerjee, 1990). 
Specifically, the following can be obtained 


nj(0<r% = U C)t£ 

(5.23a) 

*8(0 - KliO + «ft(0) = 

(5.236) 

dx j u N ,t\ _ 9L “(0 N 

d( M 

(5.23c) 


in which is obviously the nodal temperatures, and, 

Dijkl = + tySiktjl- 

Equations (5.23) form an independent set that can be resolved numerically for <rjj(£) and 
u ij(0 completely in terms of known nodal quantities u% u and t% w , without the need for kernel 
integration nor convolution. Notice, however, that shape function derivatives appear in 
(5.23c), thus constraining the representation of stress on the surface element to something 
less than full quadratic variation. The interior stress kernel functions, defined by (5.22) 
are also detailed in Appendix B. 
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6. COMPUTER PROGRAM DEVELOPMENT 

6.1 Introduction 

The goal of the computer program developed for ceramic composites is the accurate and 
efficient implementation of the formulation described in Chapters 2-5. Of equal importance 
is the degree of generality required in the definition of component geometry, loading and 
material properties. This is necessary if the program is to be applicable to real problems 
in the aerospace industry. 

For this reason the ceramic composite formulation has been implemented in the three- 
dimensional boundary element computer code ‘BEST3D’ (Boundary Element Stress Tech- 
nology - Three-dimensional) which was developed for NASA by Pratt and Whitney and 
SUNY /Buffalo under contract NAS 3-23697. Since its development, BEST3D has been 
proven to be a highly accurate and numerically efficient boundary element program. 

The development of the computer program ‘Composite-BEST’ is discussed in the fol- 
lowing sections. 

6.2 Program Structure 

Composite-BEST is designed to be a fully general ceramic composite analysis sys- 
tem employing the boundary element method. The program is written using standard 
FORTRAN 77. Development has been carried out at SUNY/Buffalo on an HP 9000 mini- 
computer, a SUN-4 workstation, and a SUN Sparcstation-1. The nature of the method 
is such that, for any realistic problem, not all required data can reside simultaneously in 
core. For this reason extensive use is made of both sequential and direct access scratch 
files. 

The program first executes an input segment. After the input has been processed, 
the surface integrals are calculated and assembled into the set of system equations using 
specified boundary conditions, followed by the insert assembly and the inclusion of the 
insert equations in the general system. The system matrix is then decomposed and saved 
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on disk, followed by the calculation of the solution vector. The full displacement and 
traction (and/or temperature and flux) solution on each boundary element and insert 
element is then reconstructed from the solution vector. In a time dependent problem the 
process of constructing the load vector for the system equations is repeated at each time 
step, but the integration, formation and decomposition of the system matrix are done 
only once. However, in an analysis with nonlinear matrix-insert interface connections, the 
assembly and decomposition must be carried out several times. 

The current program limits are as follows: 

- 20 time points (solutions for different loadings) 

- 15 generic modeling regions (GMR) 

- 600 elements (300 elements per GMR) 

- 100 inserts per generic modeling region 

- 500 inserts per problem 

- 2500 modeling nodes 

- 1200 source points (600 source points per GMR) 

Various aspects of the computer program are discussed below. 

6.3 Program Input 

The input for Composite-BEST is free field. Meaningful keywords are used to identify 
data types and to name particular data sets. The input is divided into five types: 

1. Case Control Cards 

The case control cards define global characteristics of the problem. In addition to the 
problem title, the times for multiple time steps are defined. The reading or writing of 
restart data is also defined at this point. The restart facility allows one to change 
the arrangements to fibers without recalculating the various coefficients. 
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2. Material Property Definition 

The material property input allows the definition of material properties for a variety 
of materials. The Young’s modulus can be prescribed in tabular form for a user-defined 
set of temperatures. Temperature independent values of Poisson’s ratio are also defined. 

3. Geometry Input 

Geometry input is defined one GMR (generic modeling region, or subregion) at a time. 
To initiate the input, a tag is provided to identify the GMR, a material name and reference 
temperature are defined to allow initialization of material properties. 

The next block of geometry input consists of the Cartesian coordinates of the user 
input points for the outer surface geometry definition of the composite matrix, together 
with identifiers (normally positive integers) for these geometric nodes. 

Following the definition of an initial set of nodal points, the surface connectivity of the 
outer surface of the composite matrix is defined through the input of one or more named 
surfaces. Each surface is made up of a number of elements, with each element defined 
in terms of several geometric nodes. Three sided elements, defined using six rather than 
eight geometric nodes, are used for mesh transition purposes. The terms quadrilateral 
and triangle are normally used to refer to the eight and six noded elements, although the 
real geometry represented is, in general, a nonplanar surface patch. Seven and nine noded 
elements are made available by adding a central node to the six or eight noded elements. 

Over each element the variation of displacement and traction (and/or temperature 
and flux) can be defined using either the linear or quadratic shape functions. Linear and 
quadratic elements can share a common side, which is then constrained to have linear 
displacement and traction (and/or temperature and flux) variation. 

Finally an option is available to allow quadratic functional variation (8 or 6 nodes) 
to be used in conjunction with linear geometry (4 or 3 nodes). In this case the program 
generates the additional nodes automatically at mid-point of the sides. The characteristics 
of the various element types are summarized below. 
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Geometry 


Surface Element Type Nodes 

Linear Quadrilateral 8 or 9 

Linear Triangle 6 or 7 

Quadratic Quadrilateral 8 or 9 

Quadratic Triangle 6 or 7 

Quadratic Quadrilateral 4 

Quadratic Triangle 3 


Displacement/Traction Nodes 
(and/or Temperature/Flux Nodes) 

4 

3 

8 or 9 
6 or 7 
8 
6 


Following the definition for the composite matrix outer surface, the embedded inserts 
are then defined. These are defined as curvilinear line elements with a prescribed radius 
of the cross-section. The inserts are generally straight, however as noted, curved inserts 
are also allowed. The user first defines the nodal coordinates of the centerline of the 
insert. Thereafter, the radius and the insert connectivity is defined. Linear and quadratic 
elements are available for both geometry and functional variation, however, quadratic 
functional variation over linear geometry is not presently available. The various options 
for the insert elements are summarized below. 

Geometry Displacement/Traction Nodes 
Insert Element Type Nodes (and/or Temperature/Flux Nodes) 

Linear- Linear 2 2 

Quadratic-Linear 3 2 

Quadratic-Quadratic 3 3 


Note only the surface of the insert needs to be defined, i.e., the hole in the composite 
matrix which encompasses the insert does not have to be explicitly defined. 

4. Interface Conditions 

The interface input describes the connection of surfaces or elements of one region to 
another and between the matrix and the inserts. Special types of matrix-insert interface 
conditions which are available presently include fully-bonded and sliding contact (including 
coulomb friction) and springs. 
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5. Boundary Condition Input 

The final input section provides for the definition of boundary conditions, as functions 
of both position and time. Data can be input for an entire surface, or for a subset (ele- 
ments or nodes) of a surface. Input can be in global coordinates, or can define rollers or 
pressure (or flux) in the local coordinate system. Input simplifications are available for 
the frequently occurring cases of boundary data which is constant with respect to space 
and/or time variation. Each boundary condition set can be defined at a different set of 
times. 

6.4 Surface Integral Calculation 

Following the processing of the input data, the surface integrals occurring in equations 
(2.6), (2.7), (5.6), (5.18) and (5.21) are evaluated numerically. This is the most time 
consuming portion of the analyses. In Composite-BEST the results of these integrations 
are stored as they are calculated, rather than being assembled into the final equation system 
immediately. Although this is somewhat more costly in terms of storage and CPU (central 
processing unit) time, it has led to much greater clarity in the writing of Composite-BEST. 
In addition, it provides much greater flexibility in the implementation of various restart 
and boundary condition options. 

The calculations proceed first by GMR (generic modeling region), then by source point 
(the equation being constructed) and finally by surface element and insert element. The 
results for each source point element pair are written to disk. All of the calculations are 
carried out and stored in the global (Cartesian) coordinate system. 

The integration of the BEM equations is the most complex part of the code. In this 
process either singular or nonsingular integrals can be encountered. The integrals are 
singular if the source point for the equations being constructed lies on the element being 
integrated. Otherwise, the integrals are nonsingular, although numerical evaluation is still 
difficult if the source point and the element being integrated are close together. 

In both the singular and nonsingular cases Gaussian integration is used. The basic 
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technique is developed in Banerjee and Butterfield, 1981. In the nonsingular case an 
approximate error estimate for the integral was developed based on the work of Stroud 
and Secrest (1966). This allows the determination of element subdivisions and orders of 
Gaussian integration which will retain a consistent level of error throughout the structure. 
Numerical tests have shown that the use of 3, 4 and 5 point Gauss rules provide the best 
combination of accuracy and efficiency. In the present code the 4 point rule is used for 
nonsingular integration, and error is controlled through element subdivision. The origin of 
the element subdivision is taken to be the closest point to the source point on the element 
being integrated. 

If the source point is very close to the element being integrated, the use of a uniform 
subdivision of the element can lead to excessive computing time. This frequently happens 
in the case of aerospace structures, due either to mesh transitions or to the analysis of 
thin walled structures. In order to improve efficiency, while retaining accuracy, a graded 
element subdivision was employed. Based on one- dimensional tests, it was found that 
the subelement divisions could be allowed to grow geometrically away from the origin of 
the element subdivision. Numerical tests on a complex three-dimensional problem have 
shown that a mesh expansion factor as high as 4.0 can be employed without significant 
degradation of accuracy. 

In each case of singular integration (source point on the elements being integrated) 
the element is first divided into subelements. The integration over each subelement is 
carried out using a Jacobian transformation in mapping. This coordinate transformation 
produces nonsingular behavior in all except one of the required integrals. Normal Gauss 
rules can then be employed. The remaining integral (that of the traction kernel Fij times 
the isoparametric shape function which is 1.0 at the source point) is still singular, and 
cannot be numerically evaluated with reasonable efficiency and accuracy. Its calculation is 
carried out indirectly, using the fact that the stresses due to a rigid body translation are 
zero (Lachat and Watson, 1976). It has been found that subdivision in the circumferential 
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direction of a two-dimensional surface element is required to preserve accuracy in the 
singular integration of the outer surface. A maximum included angle of 15 degrees is used. 
Subdivision in the radial direction has not been required. 

The integrals required for calculation of displacement, stress, temperature, and flux at 
interior points are of the same type as those involved in the generation of the system equa- 
tions, except that only nonsingular integrals are involved. If the source point involved is 
located on the surface of the body, then numerical integration is not required. Instead, the 
required quantities are calculated using the displacements and tractions (or temperature 
and flux) on the element (or elements) containing the source point, as discussed in Section 
2.5. 

6.5 System Matrix Assembly 

The first step in the assembly process is the reduction of the rectangular matrix of 
F integrals to a square matrix. This matrix is the prototype of the system matrix. The 
columns of the matrix are transformed or replaced, as required by the boundary conditions, 
as the assembly process proceeds. 

The next step in the process is the incorporation of the insert equations in the system. 
As was described in detail for each specific analysis in Chapters 2-5, the insert assembly 
consists of an insert by insert matrix multiplication and backsubstitution. The backsub- 
stitution minimizes the number of equations required in the system by eliminating some 
of the unknown quantities on the inserts. 

A key problem in the entire process is the proper definition of appropriate coordinate 
systems, on a nodal basis. This is a problem common in any direct boundary element 
method which treats structures with nonsmooth surfaces. It arises because the tractions 
at a point are not uniquely determined unless the normal direction to the surface varies 
continuously at the point in question. 

The original surface integral calculations are all done in global coordinates. If the 
displacement (or temperature) boundary condition is specified at a given node, in global 
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coordinates, then no new coordinate system definition is required. It is only necessary to 
keep track of the subset of elements, containing the given node, on which the fixed dis- 
placement (or temperature) is to be reacted. However, if a displacement (or temperature) 
is specified in a nonglobal direction at a given node, then a new nodal coordinate system 
must be defined and, potentially, updated as further boundary conditions are processed. 
The associated nonzero reactions must then be expressed in the new coordinate system. 

Following this preparatory work, the final assembly of the system equations is carried 
out. It is performed in three major steps: 

1. Transformation of the columns of the matrices to appropriate local coordinate systems 
and incorporation of any boundary conditions involving springs. 

2. Incorporation of compatibility and equilibrium conditions on interfaces. 

3. Application of specified displacements and tractions (and/or temperatures and fluxes). 

Two particular features of the equation assembly deserve special comment. First, 
in multi- G MR problems the system matrix is not full. Rather it can be thought of as 
consisting of an NxN array of submatrices, each of which is either fully populated or 
completely zero. Only the nonzero portions of the system equations are preserved during 
system matrix assembly. In order to improve the numerical conditioning of the system 
matrix for the solution process, the columns are reordered so that the variables from the 
two regions, lying on the same interface, axe as close together as possible. 

Second, rather than simply assembling an explicit load vector at each time point in the 
solution process, load vector coefficient matrices are assembled and stored. These allow 
the updating of the load vector at any required time point simply by interpolating the 
time dependent boundary conditions and performing a matrix multiplication. A Similar 
process is used in the calculation of interior and boundary stresses. 
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6.6 System Equation Solution 

The solver employed in Composite-BEST operates at the submatrix level, using soft- 
ware from the UNPACK package (Dongarra, 1979) to carry out all operations on subma- 
trices. The system matrix is stored, by submatrices, on a direct access file. The decom- 
position process is a Gaussian reduction to upper triangular (submatrix) form. The row 
operations required during the decomposition are stored in the space originally occupied 
by the lower triangle of the system matrix. Pivoting of rows within diagonal submatrices 
is permitted. 

The calculation of the solution vector is carried out by a separate subroutine, using the 
decomposed form of the system matrix from the direct access file. The process of repeated 
solution, required for problems with multi-time steps, is highly efficient. 

6.7 Nonlinear Solution Process 

The nonlinear solution algorithm starts with an elastic analysis (with linear elastic 
spring interface connection between the matrix and the insert) of the problem for the 
first loading increment (complete with the specified boundary and body force loading). 
At the end of the elastic increment the state variables are calculated and the nonlinear 
constitutive interface relations are established. The difference between the actual tractions 
and the elastic (or previous) tractions is estimated. A new equation system is assembled 
with the calculated nonlinear constitutive interface relations. The process is essentially 
repeated until the constitutive equations yield tractions that are negligibly different from 
the previous step. A new loading increment is taken and the process is repeated for each 
subsequent increment. 

6.8 Output Description 

The output from Composite-BEST is relatively straightforward. It consists of ten 
sections, as follows: 

1. Complete echo of the input data set. 
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2. Summary of case control and material property input. 

3. Complete definition for each GMR, including all surface insert nodes, surface and insert 
elements. 

4. Complete summary for each interface and boundary condition set, including the ele- 
ments and nodes affected. 

5. Boundary solution (on an element basis), including displacements and tractions (and/ 
or temperature and fluxes) at each node of each element. 

6. The resultant load on each element and on the entire GMR is calculated and printed. 

7. Solution for the displacements and tractions (and/or temperatures and fluxes) at the 
Insert-Matrix composite interface (on an element basis). 

8. Displacement, stress and strain (and/or temperature and fluxes) on a nodal basis, at 
all surface nodes, for each GMR. 

9. Displacements (and/or temperatures) at interior nodes. 

10. Stresses at interior nodes. 
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7. EXAMPLES OF FIBER COMPOSITE ANALYSIS 

7.1 Introduction 

In this section a number of examples are presented to verify and demonstrate the 
applications of the ceramic composite formulation of the previous chapters. 

In the mesh diagrams of the models containing the inserts, a double line is used to 
indicate the centerline of the insert elements. The length of these elements are shown in 
proper proportion for the three-dimensional views, however, the radii of the inserts are not 
indicated on these diagrams. The double line is a symbolic representation of the insert 
elements and does not in any way indicate the diameter of the insert. Refer to the example 
description for the values of the radii. 

Throughout this section consistent units are used in the definition of the examples. 
This means all lengths are defined in the same units and the tractions and the elastic 
moduli are defined in terms of these lengths as force/length 2 . No confusion should arise 
since the results are reported as non-dimensional quantities. 

7.2 Cube with a Single Insert 

The first test of the formulation is on a unit cube with a single insert of radius 0.1 
through the center of the cube. The cube is subjected to tension and shear in the direction 
parallel and perpendicular to the insert. The cube has a modulus of 100.0 and a Poisson 
ratio of 0.3. Consistent units are used for all information described in this problem. An 
insert with two different moduli of 1,000 and 10,000 is studied. The Poisson ratio of the 
insert is assumed to be the same as that of the cube. 

The problem is analyzed by both the present formulation and by a full three-dimensional 
multi-region BEM approach. As shown in Fig. 7.2.1, the model for the insert formulation 
consists of fourteen quadratic boundary elements and the insert contains three quadratic 
insert elements. The two-region, three-dimensional model shown in Fig. 7.2.2 contains 
twenty quadratic boundary elements in the first region and sixteen in the second. Note 
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9-noded elements are used in describing the insert and hole to accurately capture the 
curvilinear geometry. 

In Fig. 7.2.3, the profile of the end displacement of the cube under a uniform normal 
traction of 100.0 (in parallel with the insert) is shown. The present formulation is in good 
agreement with the full three-dimensional results for Ei/E = 10. For the case Ei/E = 100, 
the insert formulation exhibits greater stiffness than the 3-D results. This difference is 
contributed to the way the load is distributed from the insert to the composite matrix. In 
the full 3-D model, the applied traction and the resulting reactions at the fixed end act 
directly on the end of the insert. In the composite formulation, the insert is assumed not 
to intersect the boundary surface and therefore the insert is moved back slightly from the 
end of the cube. The load is therefore transferred through the composite matrix to the 
end of the insert and to its sides in a manner that is slightly different from the full 3-D 
analysis. 

In Fig. 7.2.4, the stress distribution through the center of the cube (from A to B 
as indicated in the figure) is shown. Again the results axe very good for Ei/E = 10, and 
deviates slightly from the full 3-D results in the second case. 

In Figs. 7.2.5 and 7.2.6, the lateral displacements along the side of the cube are shown 
for a cube subjected to a shear traction of 100. For the case of applied shear perpendicular 
to the insert (Fig. 7.2.5), the results for both the insert and full 3-D model show good 
agreement. Once again a slight deviation is observed for Ei/E = 100. In the case of the 
shear traction in the plane of the insert (Fig. 7.2.6) the insert has little effect on the 
displacement (as anticipated) and all results fall in close proximity. 

7.3 Lateral Behavior of a Cube with Multiple Inserts 

Existing methods of analysis of composite material based on mechanics of materials 
have been relatively successful in predicting the behavior of composite material for loading 
in the longitudinal direction. The properties perpendicular to the direction of the fibers 
are not so readily predictable by present means. The focus of the present example concerns 
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this lateral behavior. 

Four cubes (Fig. 7.3.1) with one, two, five and nine inserts are fixed with a roller 
boundary condition on one side and subjected to a uniform traction, perpendicular to the 
inserts. The material properties, given in consistent units, are 

gintert _ jqqq ^matrix _ jqq 

^.insert _ Q 3 ^matrix _ 0.3 

For the cube with one and two inserts, the boundary mesh consists of two quadratic 
surface elements on each lateral side and four elements on the top and bottom. For the 
cubes with five and nine inserts, one additional element was added to the side with the 
applied traction and to the side with the roller boundary condition. The top and bottom 
faces contain six elements to match the pattern of the sides. In all cases, each insert 
contained three one- dimensional quadratic elements. 

The profile for the end displacement of the cube with one insert and five inserts are 
shown in Figs. 7.3.2 and 7.3.3. The results are seen to be in good agreement with the 
two-dimensional results. The 2-D results are approximations since plane stress is assumed. 
The 3-D solutions for the one insert are within 2% error of the 2-D solution and within 
3% for the case of 5 inserts. 

Also shown in Fig. 7.3.4 are the average end displacements for the one, two, five 
and nine inserts. Results show good agreement with 2-D results. For one, two and five 
inserts, the solutions are within 2% error of the 2-D results and 6% for the case of nine 
inserts where the volume ratio of the inserts to the total volume is 28.2%. The result is 
also displayed in a plot of Effective Modulus vs. Insert Volume Ratio in Fig. 7.3.5. The 
effective modulus is defined as the average stress/average strain. The three-dimensional 
results follow closely to the two-dimensional solution. 
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7.4 Thick Cylinder with Circumferential Insert Supports 

The strength of a cylinder under internal pressure can be increased by adding stiff 
circumferential insert supports. In the present example, a three-dimensional, open ended 
thick cylinder are 10 and 20 respectively, the thickness is 2 and the radius of the fully- 
bonded inserts is 0.5. By using roller boundary conditions on the faces of symmetry, only 
a fifteen degree slice of the thick cylinder needs to be modeled. As shown in Fig. 7.4.1, 
sixteen eight-noded quadratic boundary elements are used to define the sides of the model, 
a nine-noded element is used on both the internal and external faces of the cylinder, and 
three insert elements are used per insert. Note, the inserts in this problem axe curvilinear 
in geometry. The elastic modulus of the cylinder is assumed to be 100, and the effect of 
inserts with five different moduli of 100, 250, 500, 750 and 1000 is studied. The Poisson 
ratio is 0.3 for both the composite matrix and insert, and the internal pressure in the 
cylinder is 100. 

Results from a multi-region, axisymmetric BEM analysis were used for comparison 
with the 3-D insert results of the present example. The axisymmetric model consists 
of twenty quadratic boundary elements on the outer surface, and six boundary elements 
per hole and per insert (Fig. 7.4.2). The radial displacement through the thick cylinder 
along the top face is shown in Fig. 7.4.3 for all five moduli. The displacement for the 
composites with low Ei/E ratios are in good agreement with the axisymmetric results, and 
diverge slightly for higher Ei/E ratios. In Fig. 7.4.4, the circumferential stress is shown 
for the same points along the top edge. This stress is smooth for the homogeneous case 
( Ei/E = l.o) and exhibits increasing fluctuations as the Ei/E ratio increases and the inserts 
take on more of the load. The circumferential stress of the 3-D insert model is in good 
agreement with the axisymmetric results for all cases. In Fig. 7.4.5, the radial stress is 
displayed for the two models. The inserts have little effect on this stress and the curves of 
the five moduli for both approaches lie in close proximity. 
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7.5 Cube with Multiple Inserts with Random Orientation 

In an attempt to analyze a material with a random fiber structure, cubes with multiple 
inserts oriented in random directions are studied. The cubes Eire of unit length and have 
four boundary elements per side (Fig. 7.5.1a). Randomly oriented fibers of variable length 
with radii of 0.05 are placed in five cubes in quantities of 5, 10, 15, 20 and 25 (Fig. 7.5.1b- 
f). Three cases of material properties are considered for each cube. The modulus of the 
composite matrix is 100 for all cases, however, the modulus of the inserts axe 500, 10,000 
and 200,000 for the three cases studied. Poisson’s ratio is uniformly 0.3 throughout. Roller 
boundary conditions are employed on three adjacent sides and a uniform normal traction 
of 100 is applied to a fourth face. 

The normal end displacement at the center of the face on the side with the applied trac- 
tion is plotted against the number of inserts in a cube for the three materials (Fig. 7.5.2). 
The displacement decreases with increasing number of inserts per cube and increasing Eij E 
values as expected. 

7.6 A Beam with Insert Reinforcement in Bending 

In the last example, the applicability of the present formulation to the study of the 
micromechanical behavior of the ceramic composite is apparent. The present formulation, 
however, is equally applicable to typical problems encountered by civil engineers. Rein- 
forced concrete can now be modeled exactly as a three-dimensional body and studied in 
detail for the first time. The present example considers a reinforced concrete beam. Here 
the concrete plays the role of the composite matrix and the reinforcement bars play the 
role of the fiber insert. In Fig. 7.6.1, a 4 x 1 x 1 beam with four inserts is modeled using 
twenty-eight quadratic boundary elements. The effect of the ratio of insert modulus to 
matrix modulus ( Ei/E ) is studied for a range of values between 1 and 100. The Poisson 
ratio is 0.3 for both the beam and reinforced rods. 

The beam is completely fixed at one end and a downward shear traction of 100 is 
applied to the other end. The non-dimensional vertical displacement of the end obtained 
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from the present analysis is shown in Fig. 7.6.2 as a function of Ei/E. The non-dimensional 
displacement is defined as the end displacement of the reinforced beam divided by the 
displacement of a homogeneous beam under similar boundary conditions. 

The end displacement obtained from the mechanics of material solution is also displayed 
in Fig. 7.6.2 in non-dimensional form. The curvature of the two plots are very similar but 
differ in magnitude. This difference is contributed to the fact that although the mechanics 
of material solution accounts for the stiffening due to the inserts, it does not include the 
effect of interaction between inserts. 

7.7 Laminated Fiber Composite 

A laminated composite fabricated from a fiber composite material is shown in Fig. 
7.7.1. The fiber composite is constructed with a single row of fully-bonded fibers oriented 
in the same direction. A two-ply laminate is then constructed from the fiber composite 
with the fibers of the two layers oriented at 90° angles. A boundary element model created 
for the study of this material is shown in Fig. 7.7.2. A small slice containing two inserts in 
each layer is used. The model consists of two regions. The outer surface of each region is 
modeled with sixteen quadratic boundary elements and each insert contains two quadratic 
insert elements. The interface between the two regions is assumed to be a perfect bond, 
however, the present version of the program also allows for sliding and spring connections. 

The composite structure is subjected to bi-axial tension. This is accomplished with 
normal tractions of 100 applied to two adjacent roller boundary conditions applied to the 
opposite ends. The elastic modulus of the composite matrix of both regions are assumed 
to be 100, and the moduli of the inserts very between 100 and 10,000. The Poisson ratio 
is 0.3 for both the composite matrix and inserts at all times. 

Figure 7.7.3 displays the displacement as a function of insert moduli for a point on the 
interface at the corner of the plate adjacent to the sides with the applied traction. The 
material exhibits less displacement as the modulus is increased, as expected. 
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7.8 Heat Conduction: Cube with Random Inserts 

The cube with randomly oriented inserts shown in Fig. 7.5.1 was analyzed for heat 
conduction analysis. The left end was subjected to a prescribed temperature of 100°F and 
for the right end a temperature of 0° was specified. All other surfaces were assumed to be 
insulated. Figure 7.8.1 shows the equivalent thermal conductivity of a cube for different 
insert arrangements. In this analysis the ratio of thermal conductivities between the insert 
and the matrix was assumed to be 100. 

7.9 Thermoelasticity: Effective Coefficient of Thermal Expansion 

The addition of fiber inserts in a material alters the thermal expansion of the material. 
The effective coefficient of thermal expansion of a composite material is dependent on many 
factors such as: The elastic and thermal properties of the individual constituents; the size, 
shape, orientation, and number of fibers; and the interface interaction between the fiber 
and matrix. Some investigators [Hopkins and Chamis, 1985] have formulated equations for 
the determination of the effective material properties of fiber composite materials, includ- 
ing the coefficient of thermal expansion. These equations axe, however, limited to specific 
types of insert arrangements and interface connections. For general arrangements, experi- 
ments can be carried out for the prediction of effective material properties. Experiments, 
however, are expensive and time consuming. The finite element method may be employed 
for this purpose, however, the convergence of the solution may be slow and the solution 
may be expensive and difficult to achieve. The present uncoupled thermoelastic BEM 
implementation is an ideal alternative for the prediction of effective material properties 
for general, tubular fiber composite materials. The analysis is relatively simple and cost 
efficient. 

A boundary element model of a cube with nine insert fibers is shown in Fig. 7.9.1. 
The fibers are assumed to be perfectly bonded to the composite matrix. Sixteen boundary 
elements are used to model the outer surface of the cube and one insert element is used to 
model each of the nine fibers. The cube is subjected to a uniform temperature increase by 
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simply specifying the temperature on one face of the cube and zero flux on the other five 
sides. The cube is allowed to expand freely, however, rigid body translation is prevented. 
The user specified radius of the insert fibers is easily changed to simulate various void 
ratio, therefore, minimizing the effort required for re-analysis of the cube with different 
insert to total volume ratios. 

In Fig. 7.9.2 the effective coefficient of thermal expansions in the lateral and transverse 
directions are shown as a function of the insert to total volume ratio for a fiber composite 
material with the following material properties: 


^matrix _ Q x lQ 6 p S i 

^insert _ 2 g Q x 10 6 pgi 

^matrix _ 0 24 

^insert _ q .24 

matrix _ q g x 10 ~6/ F 

a insert _ 3 j x IQ ~ 6 / F 

Also shown are the predictions by Hopkins and Chamis (1985). The solutions of the two 
analysis are in good agreement at low insert to total volume ratios and deviate slightly 
from one another as the insert to total volume ratio is increased. A graph for a similar 
analysis is shown in Fig. 7.9.3 for a perfectly-bonded fiber composite material with the 
following material properties of aluminum and steel are used. 

matrix _ jq q x 1Q 6 pg j 

^insert _ 3 q Q x 1Q 6 pgi 

^matrix _ q 3 

^insert _ q 3 

Q matrix _ 12 . 0 X 10 ~ 6 / F 

a insert _ g Q x IQ ~*fF 


This time the coefficient of thermal expansion of the inserts is less than that of the matrix, 
and therefore, the effective coefficient of thermal expansion decreases with increasing insert 
to total volume ratio. Once again good agreement is seen between the solutions of the two 
analyses. 

7.10 Heat Conduction: Linear Thermal Resistant Insert Interface 

A cube with five inserts is used to demonstrate the linear thermal resistant interface 
relation between the matrix and the insert. The boundary element mesh containing sixteen 
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surface elements and five insert elements is shown in Fig. 7.10.1. The radius of the inserts 
is 0.1. The conductivity of the insert and matrix are both 100.0, however, the conductivity 
across the matrix/insert interface is varied. The cube is subjected to a temperature of 
100.0 on a face parallel to the inserts and a temperature of 0.0 on the opposite face. The 
other four sides axe insulated. 

A graph of the flux at the center of the zero temperature face versus the thermal con- 
ductivity of the interface is shown in Fig. 7.10.2. For very large values of the interface 
conductivity the interface connection approaches a perfect bond and the overall solution 
approaches the solution of a homogeneous cube with a conductivity of 100.0. As the in- 
terface conductivity decreases the flux decreases. For an interface conductivity of 0.0 the 
cube exhibits the solution of a cube with perfectly insulated holes. 

7.11 Thermoelasticity: Linear Spring-Thermal Resistant Insert Interface 

The cooling of a cube containing five fiber inserts is analyzed using the thermoelastic 
formulation. Shown in Fig. 7.10.1 is a boundary element model with sixteen surface 
elements and five insert elements of radius 0.1. The inserts are assumed to be connected 
by a linear spring-thermal resistant interface. To demonstrate the effect of the interface 
connection, the insert and the matrix are assumed to be constructed from the same material 
with an elastic modulus, a Poisson ratio, a thermal conductivity and a coefficient of thermal 
expansion of 100.0, 0.3, 100.0, and 0.01, respectively. The thermal conductivity across the 
matrix-insert interface is held constant at 100.0, however, the spring constants between the 
insert and the matrix are varied. A roller boundary condition is applied on two opposite 
faces of the cube and the other four sides axe free to expand, however, rigid body translation 
is prevented. The cube is then subjected to a uniform decrease in temperature of 100.0. 

A graph is shown in Fig. 7.11.1 of the resultant tractions at the center of the roller 
boundary condition faces versus the interface spring constants ( k n = k t ). For large spring 
constant values the interface connection approaches a perfect bond, and since the material 
properties of the insert and matrix are equal, the overall solution approaches the solution of 
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a homogeneous cube. As the spring constants decrease, the tractions at the roller boundary 
conditions decrease. 

7.12 Nonlinear Matrix-Insert Interface: Beam in Bending 

A boundary element model of a 4x1x1 beam with four longitudinal inserts from example 
7.6 is shown in Fig. 7.6.1. The effect of the ratio of the insert modulus to matrix modulus 
(Ei/E) is studied for a range of values between 1.0 and 100.0. The Poisson ratio is 0.3 
throughout. A nonlinear frictional interface between the insert and matrix is assumed. 
The coefficient of friction is 0.2 and the normal tangential spring coefficients are both 1000.0 
times the elastic modulus of the matrix. 

The beam is completely fixed at one end and a downward shear traction of 100.0 is 
applied to the other end. The non-dimensional vertical displacement of the end obtained 
from the present analysis is shown in Fig. 7.12.1 as afunction of Ei/E. The non-dimensional 
displacement is defined as the end displacement of the fiber composite beam divided by 
the displacement of a homogeneous beam subjected to the same boundary conditions. The 
plot of end displacement obtained in example 7.6 for a reinforced beam perfectly bonded 
insert-matrix interface connections is also displayed in Fig. 7.12.1. The beam with the 
perfectly-bonded inserts exhibits greater stiffness than the beam with the nonlinear insert 
interface connections, as expected. 

7.13 Effect of Poisson Ratio in Fiber Composites 

The goal of the present five year project is to develop an accurate, yet practical bound- 
ary element program for the analysis of ceramic composite material. A conventional BEM 
approach is computationally expensive and mesh generation is tedious for an analysis of 
a solid which may have hundreds of fiber inserts. Therefore the challenge of the present 
work is to introduce certain approximations in the formulation which will produce an effi- 
cient BEM analysis and will not compromise accuracy. One approximation unique to this 
BEM implementation is the assumption that boundary quantities vary as a trigonomet- 
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ric function around the circumference of the insert and hole. This was shown to be an 
accurate approximation through the numerical applications of this section. Secondly, the 
steady-state BEM formulations derived in Chapters 2-4 assume that the Poisson ratio of 
the insert is equal to the Poisson ratio of the matrix. (This limitation was removed in the 
formulations of Chapters 5, however, these formulations are more expensive.) This second 
assumption is investigated in this example using a two-dimensional BEM analysis. 

A boundary element model of a cube with a single insert is shown in Fig. 7.13.1. 
Sixteen quadratic boundary elements are used to model the outer boundary and eight 
elements are used to model both the hole and insert. The modulus of elasticity of the 
matrix is 10 6 and the Poisson ratio is 0.25. The elastic modulus of the insert is assumed 
to range from 10 4 to 10 s and the Poisson ratio from 0.0 to 0.5. The cube is secured by 
roller boundary conditions on two adjacent sides, a traction of 10 5 is applied on a third 
side, and the normal displacement is measured at the center of the remaining face. The 
ratio of insert diameter to side of cube is 0.2 for the volume of insert to total volume ratio 
of 0.031. The displacement versus the variation in Poisson ratio (^naerty^matrix^ j s shown in 
Fig. 7.13.2 for the elastic modulus ratios (£ in * ert /^matrix j c f o.Ol, 0.1, 1.0, 10.0, and 100.0. 
The displacement in these curves are normalized by the displacement obtained for the 
respective elastic modulus ratios when the Poisson ratio of the insert and matrix are the 
same (v 1 /v M = 1.0). The displacement is observed to be minimal in all cases. The effect is 
greatest for the modulus ratio E l jE M of 1.0 at the extreme values of u 1 /v M . The maximum 
error, however, is less than 4.5%. 

This same problem is analyzed for the insert diameter to side of cube ratio of 0.6 for 
an insert volume to total volume ratio of 0.283 as shown in Fig. 7.13.3. The normalized 
displacement versus v 1 jv M for five modulus ratios is shown in Fig. 7.13.4. Very minute 
change with variation of v l jv M is observed for modulus ratios of 0.01 and 100.0 which are 
almost indistinguishable from one another on the graph. The variation in Poisson ratio 
is seen to have a modest effect on the displacement values for modulus ratios of 0.1 and 
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10.0. For an elastic modulus of 1.0, the effect of Poisson ratio is somewhat tolerable for 
the midrange values of i s 1 /v M , however, it has a strong effect for extreme values of v 1 fv M . 

Next the effect of the Poisson ratio in a cube with multiple fibers interacting with one 
another is investigated. A cube with nine inserts (Fig. 7.13.5) is analyzed under similar 
conditions as the analyses of a cube with a single fiber. The volume of insert to total 
volume ratio is again 0.283. In Fig. 7.13.6 the displacement versus u 1 jv M is shown for the 
multiple fiber composite for a modulus ratio E I /E M of 1.0. For comparison, the curve for a 
cube with a single insert with the same modulus and insert to total volume ratio is shown. 
The effect with the multiple fiber inserts is just slightly less than the effect with a single 
insert. 

In conclusion, we note that a difference in Poisson ratio between an insert and a matrix 
should not be ignored if the difference in the modulus of elasticity is small and the volume of 
insert to total volume is large. Nevertheless, fiber inserts are generally added to material in 
order to strengthen the material, and the elasticity modulus ratio (E 1 /E M ) is usually greater 
than five. Furthermore, many composite constituents often have comparable Poisson ratio 
values. Therefore, the assumption that the Poisson ratio of the insert must be equal to 
the Poisson ratio of the matrix is not very limiting. 

7.14 Composite with Fifty-one Fiber Inserts 

In Fig. 7.14.1 a boundary element mesh is shown for the outer surface of a rectangular 
block containing fifty small insert fibers and one large fiber (insert elements not shown). 
The insert fibers are aligned parallel to one another in the arrangement shown in Fig. 
7.14.2. Due to the relative size of the large central fiber, this fiber is explicitly modeled as 
a separate region in the conventional boundary element sense. The other fifty fibers are 
each modeled with one quadratic insert element. The total number of nodes in the problem 
is five-hundred and ninety-eight. The material properties of the composite constituents 
are 
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Matrix 


Fiber Type 1 Fiber Type 2 


E M = 17.2 x 10 6 psi 

v M = 0.24 

a M = 0.5 x 10 6 /F 


E n = 60.0 x 10 6 psi 

u 11 = 0.24 

Q n = 4.0 x 10- 6 /F 


E n = 28.0 x 10 6 psi 

v n = 0.24 

a /2 = 2.1 x 10- 6 /F 


In Fig. 7.14.3 through 7.14.9 displacement contours are shown for a model constrained 
with roller boundary conditions on three adjacent sides (bottom, left, and back faces) and 
subjected to a normal traction of 10 5 on the small (right side) face perpendicular to the 
fibers. In. Figs. 7.14.4, 7.14.6, 7.14.8 and 7.14.9 a uniform temperature decrease of 2400°F 
has also been applied in addition to the normal traction. Contours for displacement (of 
the top, traction free face) in the direction parallel to the normal traction is shown in 
Figs. 7.14.3 and 7.14.4. In Fig. 7.14.5 and Fig. 7.14.6 the contours are shown for the 
displacement (of the top face) in the direction perpendicular to the normal traction and 
the fibers. In Fig. 7.14.7 and Fig. 7.14.8 the contours are shown for the displacement (of 
the top face) parallel to the fibers. Fig. 7.14.9 displays a three-dimensional view of 7.14.8. 

In Fig. 7.14.10 a contour plot of displacement is shown for a model subjected to a 
normal traction on the top face perpendicular to the fibers. The displacement shown is of 
the top face in the direction parallel to the loading. The opposite face is fixed on rollers, 
and the lateral faces are traction free (rigid body translation is prevented). The fifty small 
insert fibers primarily strengthen the matrix uniformly, and the large insert interacts with 
the strengthened matrix resulting in the deformation displayed in the figures. 

7.15 Transient Thermoelastic Analysis of a Cube with Inserts 

An uncoupled thermoelastic analysis of a cube with one and nine inserts is studied. A 
unit cube is discretized with six quadratic boundary elements on the surface of the cube 
and one insert element per insert. The radius of the insert is 0.2 for the cube with a single 
insert and 0.1 for the cube with nine inserts. The material properties, given in consistent 
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units axe 



^insert _ jqq q 

^•matrix ^ q 


^insert _ Q.3 

^matrix _ q 3 


jj. insert _ 1000 .0 

jj, matrix _ 1Q 


pc jnsert _ 1Q 

^matrix _ 10 


The temperature of the cube is initially at zero degrees. On a face parallel to the length of 
the insert, the temperature is held at 0° degrees and the face is supported by rollers. On 
the opposite face the temperature is suddenly raised to 100°. All other sides axe insulated 
and rigid body translation is prevented. The results of the flux and end displacement 
versus nondimensionalized time for the traction free face is given in Fig. 7.15.1 through 
7.15.4. 

The three-dimensional fiber composite BEM results are compared with conventional, 
two-dimensional BEM results. In Figs. 7.15.1 and 7.15.2 the flux on the free face is 
shown for a cube with a single insert and a cube with nine inserts, respectively. Similarly, 
Figs. 7.15.3 and 7.15.4 display the results for end displacement versus time. The fiber 
composite analysis is in good agreement with the conventional, two-dimensional results. 
Slight deviation is expected due to the 2-D assumption and the use of a somewhat crude 
3-D discretization. 

7.16 Transient Heat Conduction Analysis of a Turbine Blade 

A boundary element discretization of a turbine blade is shown in Fig. 7.16.1. Half 
symmetry is employed in this model which consists of a single region with ninety-two 
quadratic boundary elements. The model is 58.2 mm long, 13.9 mm wide, the radius of 
the base is 6.95 mm, and the tip is 1.98 mm (from the axis of symmetry) in thickness at 
the largest point. 

A transient heat conduction analysis is first carried out on a homogeneous blade with 
a conductivity of ib = 0.0216. The blade initially rests in thermodynamic equilibrium at zero 
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temperature. Then, a gas at a temperature of 1200° is assumed to flow over the blade 
while a gas at a temperature of 500° surrounds the base of the blade. At the leading edge 
of the blade the film coefficient is h = 0.003395 and tapers off to h = 0.00064 at the trailing 
edge. At the base of the blade the film coefficient of h = 0.00005 is assumed. Figures 7.16.3a 
through 7.16.7a contain contour plots describing the temperature flow through the blade 
at 1.0, 2.0, 4.0, 6.0, and 8.0 seconds. 

The blade is reanalyzed, this time with eight fibers (four fibers per side) running from 
the tip of the blade to half way through the base. The radius of the fibers is 0.15 mm and 
their cross-sectional locations are shown in Fig. 7.16.2. Five (9-noded) insert elements are 
used to model each fiber. The conductivity of the fibers is 100 times the conductivity of 
the blade. In Figs. 7.16.3b through 7.16.7b the resulting temperature distributions are 
shown in contour plots. Overall, the higher conductivity of the fibers increases the heat 
flow through the blade. Hence, the heat at the tip of the blade is carried to the base faster 
than in the homogeneous blade. This results in a temperature distribution which is lower 
than the homogeneous blade near the tip and higher near the base. 
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Fig. 7.2.1 Discretization of an Insert in a Unit Cube Utilizing 
Quadratic Insert Elements 



Fig. 7.2.2 Full Three-dimensional, Multi-region Discretization 
of an Insert in a Unit Cube 
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Fig. 7.2.3 Comparison of Displacement Profiles Between the Full 
3-D Model and the Insert Element Model for the End of 
a Cube in Tension 



Fig. 7.2.4 


Axial Stress Through the Cross Section of a Unit Cube 
in Tension with a Single Insert in Parallel with the 
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Fig. 7.3.1 Arrangement of Multiple Inserts in a Unit Cube 
Subjected to Lateral Tension 
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X 



Fig. 7.3.4 Average End Displacement of a Cube Under Tension 
vs. the Volume of Insert to Total Volume Ratio 



Fig. 7.3.5 Effective Transverse Modulus of a Cube as a 
Function of Insert Volume to Total Volume 
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Fig. 7.4.3 Radial Displacement Through a Pressurized Thick 

Cylinder with Circumferential Inserts for Ei/E = l.C 
2.5, 5.0, 7.5, 10.0 



Fig. 7 .4.4 Circumferential Stress Through a Pressurized Thick 

Cylinder with Circumferential Inserts for Ei/E =1.0 
2.5, 5.0, 7.5, 10.0 
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EXFCT MODELING OF INSERTS CRXI5YM-ETRIC BEM) 
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Fig. 7.4.5 Radial Stress Through a Pressurized Thick Cylinder 
with Circumferential Inserts for Ei/E = 1 0, 2 5, 
5.0, 7.5, 10.0 
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END DISPLACEMENT / LENGTH 



NUMBER OF INSERTS 


Fig. 7.5.2 End Displacement of a Unit Cube With Random Oriented 
Inserts of Ei/E = 5.0, 10.0 and 100.0 



Fig. 7.6.1 Discretization of a Reinforced Beam Utilizing 

Quadratic Insert Elements to Model the Four Inserts 
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END DISP. /HOMOGENEOUS END DISP. 



Ei / E 

Fig. 7.6.2 Non-dimensionalized Vertical End Displacement of a 
Reinforced Beam in Bending vs. the Modulus of the 
Insert Over Modulus of the Beam 
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Fig. 7.9.2 Effective Coefficient of Thermal Expansion 

in a Cube with Nine Perfectly-bonded Insert Fibers 
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Fig. 7.9.3 Effective Coefficient of Thermal Expansion 
in a Model of an Aluminum Cube with Nine 
Perfectly-bonded Steel Inserts 
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Fig. 7 .10.1. Discretization of a Cube with Five Inserts 
Connected to the Matrix Through a Linear 
Constitutive Model 


° BErt COMPOSITE SOLUTION 


T- loo O O T - 

o * 

o o 


MRTR IX/ INSERT INTERFRCE CONDUCTIVITY 


Fig. 7.10.2 Effect of a Thermal Resistant Matrix/Insert 
Interface on the Flux Through a Fiber 
Composite 
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Fig. 7.11.1 Effect of a Spring Interface Between the 

Matrix and Insert in a Cube Subjected to a 
Uniform Temperature Decrease 



Ei ✓ E 


Fig. 7.12.1 Vertical End Displacement of a Beam in 
Bending with a Nonlinear, Frictional 
Interface Between Four Insert Fibers and 
the Matrix 
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Fig. 7.13.1. Two-dimensional, Multi-region BEM Model of 
a Cube with a Single Insert. 

(Insert Volume/Total Volume = 0.031) 
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Fig. 7.13.2. Effect of Poisson Ratio on the Lateral 

Displacement of a Single-fiber Composite in 
Tension for Values of eVe m of 0.01, 0.1, 
1.0, 10.0, and 100.0 




Fig. 7.13.3. Two-dimensional, Multi-region BEM Model of 
a Cube with a Single Insert 
(Insert Volume/Total Volume = 0.283) 
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Fig. 7.13.4. 


Effect of Poisson Ratio on the Lateral 
Displacement of a Single-fiber Composite in 
Tension for values of El/E M of 0.01, 0.1, 1.0 
10.0, and 100.0 




NORMALIZED VERTICAL DISPLACEMENT 


Fig. 7.13.5. Two-dimensional, Ten Region BEM Model of a 
Cube with Nine Inserts 
(Insert Volume/Total Volume = 0.283) 
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Fig. 7.13.6. Effect of Poisson Ratio on the Lateral 
Displacement of a Single-fiber Versus a 
Multi-fiber Composite in Tension with 

eVe m = l.o 




Fig. 7.14.1 Discretization of the Outer Surface of a 

Two-region Model Containing 50 Small Insert Fibers 
(Small Insert Fibers not shown) 


Fiber Type 
2 



Fig. 7.14.2 Cross-section of the Composite Displaying 
the Arrangement of One Large (Type 1) Fiber 
and Fifty Small (Type 2) Fibers 






Fig. 7.14.3. Contours of Displacement in the Horizontal 
Direction Relative to the Page Due to a 
(Horizontal) Traction Loading 



Fig. 7.14.4 Contours of Displacement in the Horizontal 
Direction Relative to the Page Due to a 
(Horizontal) Traction and Temperature Loading 
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Fig. 7.14.5. Contours of Displacement in the Vertical 
Direction Relative to the Page Due to a 
(Horizontal) Traction Loading 



Fig. 7 .14.6. Contours of Displacement in the Vertical 
Direction Relative to the Page Due to a 
(Horizontal) Traction and Temperature Loading 
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Fig. 7.14.7. Contours of Displacement in the Out of Plane 

Direction Due to a (Horizontal) Traction Loading 



Fig. 7,14.8. Contours of Displacement in the Out of Plane 
Direction Due to a (Horizontal) Traction and 
Temperature Loading 
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Fig. 7.14.9 Three-dimensional View: Contours of Displacement in 

the Out-of-Plane (z) Direction Due to a Traction 
(x-direction) and Temperature Loading 



Fig. 7.14.10. Contours of Displacement in the Out-of-Plane 

Direction Due to a Traction in the Out-of-Plane 
Direction (Symmetric Boundary Conditions) 
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Fig. 7.15.1 Flux vs. Time: Transient Thermoelastic Analysis 

of A Cube with One Insert 



Fig. 7.15.2 End Displacement vs. Time: Transient 

Thermoelastic Analysis of a Cube with 
One Insert 





TIME 


Fig. 7.15.3 Flux vs. Time: Transient Thermoelastic Analysis 

of a Cube with Nine Inserts 



TIME 


Fig. 7.15.4 End Displacement vs. Time: Transient Thermoelastic 

Analysis of a Cube with Nine Inserts 
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Fig. 7.16.1 Discretization of a Turbine Blade 
(Insert Elements Not Shown) 



Fig. 7.16.2 Cross-Sectional View of the Turbine Blade 

Displaying the Location of the Eight Fibers 
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Fig. 7.16.7 Transient Heat Conduction Analysis of a Turbine 
Blade at 8 Seconds 
(Contour Plots of Temperature) 



8. CURRENT ACHIEVEMENT 

8.1 Third Year Development 

The main focus of the work of this past year was the implementation of the transient 
heat conduction analysis and the implementation of the transient thermoelastic analysis. 
These analyses are vastly different from their steady-state counterparts. This and other 
important developments of the past year are outlined below. 

New Numerical Integration for Inserts : The complexity of the transient kernel func- 

tions prohibits the analytic integration about the circumference of the insert in a manner 
similar to the steady-state case. Instead a numerical integration scheme had to be devel- 
oped for this purpose. The circular shape function which was developed during the early 
stage of this work is still employed in describing the field quantities across the boundary of 
the inserts and the holes encompassing the inserts. Since numerical integration adds to the 
cost of the analysis and is required at each step in a transient BEM analysis, an efficient 
numerical integration had to be developed. The present scheme employs one-dimensional 
Gaussian integration formulas independently in two directions. This allows the order of 
the formulas to vary independently in the two directions and the accuracy and the cost 
can be optimized. 

New Efficient Assembly for Transient Analysis : The modified boundary integral ap- 

proach and the resultant assembly scheme which are used in the steady-state composite 
analyses can not be used in the transient case, since it requires the transient material prop- 
erties of the fibers to be identical to the properties of the composite matrix. For general 
fiber material properties a new approach had to be developed. The success of the steady- 
state composite analyses is based on the fact that the size of the system equations is kept 
to a minimum. This allows large problems with many inserts to run efficiently. The new 
assembly scheme for transient analysis employs backsubstitution of the insert equations 
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in an efficient manner so that the final system is the same size as the steady-state case. 
Furthermore, this backsubstitution, as well as the decomposition of the system matrix is 
only required in the initial time step. 

Alternative Method for Steady-State Comp osite Analysis: The assembly developed for 
the transient composite analysis was extended back to the steady-state formulations. The 
object of this is two-fold. First, since the assembly schemes are completely different, it 
provides a check of both methods. Second and most importantly, the original modified 
boundary integral approach is based on an assumption that the Poisson s ratio of the fibers 
are the same as the composite matrix. Although it has been shown (example 7.13) to be 
a good approximation, the new assembly scheme permits the Poisson’s ratio of the fibers 
and the composite matrix to be different. This new approach, however, is more expensive 
than the original method. 

New Nonlinear Interface Model : The first nonlinear fiber to matrix interface model 

used in the program was based on the Coulumb’s law of friction. An interface based on 
this model, however, cannot withstand any tension. If subjected to tension, a gap opens 
up and the tractions on the interface become zero. In most applications, the fiber has 
an adhesion to the composite matrix. Therefore, an interface subjected to tension will 
exhibit elastic and/or nonlinear behavior before a failure takes place. A modified version 
of the spring- Coulumb interface model has been added to the program. This model has 
the ability to sustain a traction force in both the normal and tangential directions up to a 
user specified limit. Above this limit the normal spring- Coulumb model prevails. 

Maintenance and Testing : Considerable effort has been devoted to maintaining the 

computer program and improving the quality and reliability of the code. Verification 
problems from previous years are regularly executed to assure their continued operation. 
Furthermore, the present work was thoroughly tested and verified. 
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8.2 Summary 

Significant progress has been achieved towards the goal of developing a general purpose 
boundary element program for the micro- and macro- mechanical studies of advanced ce- 
ramic composites. In the first year of the contract a ceramic composite code was developed 
based on the advanced boundary element program BEST3D. At that time the formula- 
tion for elastostatic analysis of fully-bonded inserts was implemented in the code. All the 
general purpose features in BEST3D, such as the definition of local boundary conditions 
and multi-region substructuring, were retained in order to facilitate the analysis of real 
problems encountered in industry today. 

The ceramic composite formulations for steady-state heat conduction and steady-state 
uncoupled thermoelasticity were developed and implemented in the elastostatic code de- 
veloped in the first year. Analytic integration of the kernel functions was performed about 
the circumference of the insert, and therefore, the numerical integration of the insert (and 
hole) was reduced to an efficient, one- dimensional integration along the length. An ef- 
ficient assembly and solution algorithm was developed for these formulations similar to 
the algorithm used in the elastostatic analysis. This insert assembly scheme significantly 
reduced the number of unknowns in the system, and therefore, it reduced the number of 
equations and the work required to solve them. Also, the displacement and the heat con- 
duction equations of the thermoelastic analysis are integrated simultaneously, but solved 
separately for efficiency. 

Considerable effort has been focused on the development of nonlinear interface connec- 
tions between the fibers and the composite matrix. Up to this point only perfectly-bonded 
connections were considered. The ceramic composite formulations were reformulated to 
allow the implementation of a variety of sophisticated interface relations such as spring 
connections, sliding connections with smooth or frictional resistance, progressive debond- 
ing, yielding on interfaces, and temperature dependent connections. In addition to the 
fully-bonded insert interface, constitutive relations were derived for spring connections 
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and slide connections with spring-coulomb frictional resistance. These models were im- 
plemented in a general manner to facilitate the addition of other interface models to be 
derived in the future. Since the interface relations are nonlinear and path dependent, an 
efficient nonlinear algorithm had to be derived. The nonlinear algorithm was implemented 
in both the elastostatic and the uncoupled thermoelastic analyses. 

Overall the cost of the steady-state ceramic composite analysis is just slightly more 
expensive (for a moderate number of inserts) than the cost of analyzing the matrix without 
the inserts present. The additional cost is primarily attributed to the integration of the 
outer surface of the matrix for the additional nodes on the surface of the hole containing the 
insert, and towards the expense of solving a slightly larger system. More importantly, the 
present analysis is fax more efficient than the conventional three-dimensional, multi-region 
approach which requires significantly more nodes in the discretization of the insert and hole. 
The additional nodes require additional equations which add to the expense of integration 
and solution. Furthermore, the conventional approach requires a two-dimensional numer- 
ical integration over the surface of the insert and hole as opposed to a one-dimensional 
integration used in the present method. The data preparation for the present method is 
less involved than the conventional multi-region approach and the location of the inserts 
mu be readily changed in a re-analysis. Moreover, the code developed in the present work 
allows up to 500 (100 per GMR) insert elements in an analysis. An ordinary multi-region 
code would require a 500 region capability in order to compete. In terms of computer 
expense and the cost of data preparation for a 500 region problem, such analyses would 
be impractical. 

The implementation of the transient heat conduction and the transient thermoelastic 
composite analyses are more complicated than of the steady-state cases due to the presence 
of the convolution integrals. Since the transient kernels cannot be integrated analytically, a 
n um erical integration scheme was developed for the integration of the inserts. Furthermore 
a completely new assembly scheme had to be developed for the analysis. This assembly 
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scheme, will also be adopted for the elastodynamic formulations. Although, the transient 
algorithm costs more than its steady-state counterpart, the present implementation is very 
efficient relative to the complex nature of the problem. 

The ceramic composite formulation for elastostatics, steady-state and transient heat 
conduction, and steady-state and transient uncoupled thermoelasticity has been success- 
fully implemented and has been proven to be accurate and efficient. Furthermore, the 
development and implementation of the nonlinear interface constitutive models was suc- 
cessfully carried out in a general manner. Since the boundary element method has already 
been proven successful in elastodynamic, free vibration, and failure analyses, coupled with 
the success of the present work, the plan to extend the ceramic composite formulations to 
these other areas holds great potential. In fact, the boundary element method may not 
only be the best tool, but may also be the only practical tool for the analysis of large 
ceramic composite problems encountered in industry. 
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9. FUTURE DEVELOPMENT 

Significant progress has been achieved during the first three years in developing a 
user friendly, ceramic composite program. The infrastructure of the program has been 
set up to facilitate future development. Presently the elastostatic, the steady-state heat 
conduction, and the steady-state uncoupled thermoelastic formulations for the analysis of 
fiber composite bodies with linear and frictional fiber-matrix interface connections have 
been implemented and tested. Likewise, the transient heat conduction and thermoelastic 
analysis, with bonded interfaces have been implemented and are working successfully. 
Delivery of the present code is planned for Spring of 1991. 

The remaining work for the last two years of this contract is concentrated in two areas: 
the yielding and failure of the fiber composites under loading; and the dynamic analysis 
of fiber composites. The upcoming year will focus on the first topic. An algorithm for 
yielding and tearing of individual fibers will be introduced in the program and extended 
to the yielding and failure of the composite matrix. Both a plastic fracture and a creep 
model will be included. To complete the failure analysis the nonlinear fiber-matrix interface 
development will be added to this algorithm. The end result will allow the analysis of very 
sophisticated modes of failure. 

Finally, a time saving feature will be added to the insert integration schemes which 
will eliminate some unnecessary integration when there are two or more fibers having the 
same material properties, geometry, and discretization. 

In the fifth and final year of the contract a free vibration analysis for fiber composites 
will be implemented along with a formulation for the steady-state and transient elasto- 
dynamic analysis. Next, the nonlinear interface, creep and failure models will be added 
to these implementations. Finally, the last piece of work will be the development of a 
temperature dependent interface and creep model (thermoplastic creep model) that can 
be used in conjunction with the uncoupled thermal formulation. The order of work for the 
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last two years of the contract is summarized in Table 9.1. 

The final version, in particular, will provide a very precise, yet very efficient, use- 
friendly, design and analysis tool for ceramic composites exposed to severe operating 
environments. The resulting computer code will enable an engineer to undertake rapid 
numerical experiments to gain insight into the micro- and macro-mechanical behavior of 
ceramic composite materials. Armed with this information, the disposition of fibers can 
be selected to optimize performance under inelastic, thermal and dynamic loading. 
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TABLE 9.1 




Planned Development for the 
Boundary Element Ceramic Computer Code 


Fourth Year (1991) 

- Implement a time saving insert integration feature. 

- Im plement nonlinear interfaces in transient problems 

- Implement yielding and failure of fibers 

- Implement yielding and failure of the composite matrix (plasticity fracture models, 
creep model) 

- Combine nonlinear interfaces with failure analysis 

Fifth Year (1992) 

- Implementation of steady-state elasto dynamics 

- Implementation of transient elastodynamics 

- Implementation of free vibration analysis 

- Add nonlinear interfaces and failure models to dynamic problems 

- Add temperature dependent' interfaces and failure models (thermoplastic creep model) 
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APPENDIX A - STEADY-STATE KERNEL FUNCTIONS 

This appendix contains the three-dimensional, steady-state kernel functions utilized in 
Chapters 2-4. 

A.l Heat Conduction 


For the temperature kernel, 


whereas, for the flux kernel, 
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A. 2 Elastostatics 


For the displacement kernel, 
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whereas, for the traction kernel, 
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A.3 Steady-State Thermoelasticity 


For the displacement kernel, 
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APPENDIX B - TRANSIENT KERNEL FUNCTIONS 

This appendix contains the detailed presentations of the three-dimensional, transient 
kernel functions utilized in Chapter 5. For the time-dependent problems, the following 
relationships must be used to determine the proper form of the functions required in the 
boundary element discretization. That is, 

GtyX - 0 = G oP (X - C nAt) for n = 1 

G n ap {X - 0 = G aP {X - C nAt) - G aP {X ~Z,(n- 1)A f) for n > 1, 

with similar expressions holding for all the remaining kernels. In the specification of these 
kernels below, the arguments (X ~(,t) are assumed. 

Note that the indices 

i,j, k, l vary from 1 to 3 

a,p vary from 1 to 4 
6 equals 4 

Additionally, 

Xi coordinates of integration point 

& coordinates of field point 
y« = xi-£i r 2 = yiyi. 

For the displacement kernel, 

C « = s(i(I7!rt)[(r)' <W ] 
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whereas, for the traction kernel 
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For the interior stress boundary kernels, 


where 


_ 2 (iv dFpi 


Dpij = 




+ ( 


(9G pi 

dG Pj \ 

V 9(j + 

d(i ) 


dF Pj \ 

V <90 + 

dti ) 

/ ZyiVkyk 

1 

1 


16?rr 2 ^(1 

+ 


r ; 


dG ej J_ 

< 9 £jt 87 rr V<:(A + 2 /i) 


( 3 _ 4l/ ) 

' )[(*?-%»)*«-(*?)*] 


122 
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and the prime, represents a derivative with respect to tj. Thus, 
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APPENDIX C - LIST OF SYMBOLS 


A short list of notation used in this report is given below. Bold print is used to identify 
vectors and matrices. 

A b , B b Boundary system matrices in assembled form 
Cij(() A tensor dependent on location of the field point £ 

di Difference in displacement (and/or temperature) across the matrix-insert interface 
Gij, Fij Kernel functions 

k^ Nonlinear interface constitutive relationship 
L p {m,V 2 ) Two-dimensional shape function 
N One-dimensional shape function 
M a (6) Circular shape function 
t Time 
ti Traction 
Ui Displacement 

x,X Refers to global coordinates of an integration point 
x System vector of unknown boundary quantities 
y System vector of known boundary quantities 
Cij Strain 

t] Refers to local coordinates of an integration point 
£ refers to coordinates of a field point 
<7ij Stress 

Subscript 

i,j,k Indicial notation 

i, j, k = 1 for heat conduction analysis (Chapter 2) 
i, j, k = 1,2,3 for elastostatic analysis (Chapter 2) 
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i,j,k = 1,2, 3, 4 for thermoelastic analysis (Chapters 3 and 4) 

<*,/? = 4 for heat conduction (Chapter 5) 
a, /3 = 1 , 2 , 3, 4 for thermoelasticity (Chapter 5) 

Superscript 

H refers to quantities on the surface of the hole (left for insert) in the ceramic matrix 
1, 12 refers to quantities on the surface of the (insert) fiber 

O refers to quantities on the outer surface of the ceramic matrix 
. (Incremental) time rate of change 
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